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

update scripts for adiab ions

parent 000bf3ea
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1110 additions and 313 deletions
matlab/assemble_cosolver_matrices.m
+ 30
− 30
View file @ 79443f28
codir = '/home/ahoffman/cosolver/'; codir = '/home/ahoffman/cosolver/';
matname= 'gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5'; matname= 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8';
matdir = [codir,'SGGK_tau1e-2_P4_J2_dk_0.05_kpm_5.0/matrices/']; matdir = [codir,'LDGKii_P16_J9_NFLR_8/matrices/'];
kperp = load([matdir,'kperp.in']); kperp = load([matdir,'kperp.in']);
Nkp = numel(kperp); Nkp = numel(kperp);
outdir = '/home/ahoffman/gyacomo/iCa/'; outdir = '/home/ahoffman/gyacomo/iCa/';
...@@ -13,34 +13,34 @@ end ...@@ -13,34 +13,34 @@ end
Nmax = numel(kperp); Nmax = numel(kperp);
for n_=1:Nmax for n_=1:Nmax
n_string = sprintf('%5.5d',n_); disp(n_string); n_string = sprintf('%5.5d',n_); disp(n_string);
filename = ['ei.',n_string,'.h5']; filename = ['self.',n_string,'.h5'];
if(exist([matdir,filename])>0) if(exist([matdir,filename])>0)
% Load matrices and write them in one h5 file % Load matrices and write them in one h5 file
olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5']; % olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
load_write_mat(infile,olddname,outfile,newdname); % load_write_mat(infile,olddname,outfile,newdname);
%
olddname = '/Ceipj/CeipjT'; infile = [matdir,'ei.',n_string,'.h5']; % olddname = '/Ceipj/CeipjT'; infile = [matdir,'ei.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
load_write_mat(infile,olddname,outfile,newdname); % load_write_mat(infile,olddname,outfile,newdname);
%
olddname = '/Ciepj/CiepjT'; infile = [matdir,'ie.',n_string,'.h5']; % olddname = '/Ciepj/CiepjT'; infile = [matdir,'ie.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
load_write_mat(infile,olddname,outfile,newdname); % load_write_mat(infile,olddname,outfile,newdname);
%
olddname = '/Ciepj/CiepjF'; infile = [matdir,'ie.',n_string,'.h5']; % olddname = '/Ciepj/CiepjF'; infile = [matdir,'ie.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
load_write_mat(infile,olddname,outfile,newdname); % load_write_mat(infile,olddname,outfile,newdname);
%
olddname = '/Caapj/Ceepj'; infile = [matdir,'self.',n_string,'.h5']; % olddname = '/Caapj/Ceepj'; infile = [matdir,'self.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
mat_ee = load_write_mat(infile,olddname,outfile,newdname); % mat_ee = load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5']; olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname]; newdname = ['/',sprintf('%5.5d',n_-1),olddname];
mat_ii = load_write_mat(infile,olddname,outfile,newdname); mat_ii = load_write_mat(infile,olddname,outfile,newdname);
mat_ee = 0;
% to verify symmetry % to verify symmetry
verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii)))); 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; if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
...@@ -51,16 +51,16 @@ for n_=1:Nmax ...@@ -51,16 +51,16 @@ for n_=1:Nmax
break break
end end
end end
olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5']; olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
Pmaxe = h5readatt(infile,olddname,'Pmaxe'); % Pmaxe = h5readatt(infile,olddname,'Pmaxe');
Jmaxe = h5readatt(infile,olddname,'Jmaxe'); % Jmaxe = h5readatt(infile,olddname,'Jmaxe');
Pmaxi = h5readatt(infile,olddname,'Pmaxi'); Pmaxi = h5readatt(infile,olddname,'Pmaxi');
Jmaxi = h5readatt(infile,olddname,'Jmaxi'); Jmaxi = h5readatt(infile,olddname,'Jmaxi');
dims_e = [0 0]; % dims_e = [0 0];
dims_e(1)= Pmaxe; dims_e(2) = Jmaxe; % dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
dims_i = [Pmaxi; Jmaxi]; dims_i = [Pmaxi; Jmaxi];
h5create(outfile,'/dims_e',numel(dims_e)); % h5create(outfile,'/dims_e',numel(dims_e));
h5write (outfile,'/dims_e',dims_e); % h5write (outfile,'/dims_e',dims_e);
h5create(outfile,'/dims_i',numel(dims_i)); h5create(outfile,'/dims_i',numel(dims_i));
h5write (outfile,'/dims_i',dims_i); h5write (outfile,'/dims_i',dims_i);
% %
......
matlab/load/load_params.m
+ 11
− 3
View file @ 79443f28
function [DATA] = load_params(DATA,filename) function [DATA] = load_params(DATA,filename)
DATA.inputs.CO = h5readatt(filename,'/data/input/coll','CO'); DATA.inputs.CO = h5readatt(filename,'/data/input/coll','CO');
DATA.inputs.K_N = h5readatt(filename,'/data/input/ions','k_N');
DATA.inputs.K_T = h5readatt(filename,'/data/input/ions','k_T');
DATA.inputs.Q0 = h5readatt(filename,'/data/input/geometry','q0'); DATA.inputs.Q0 = h5readatt(filename,'/data/input/geometry','q0');
DATA.inputs.EPS = h5readatt(filename,'/data/input/geometry','eps'); DATA.inputs.EPS = h5readatt(filename,'/data/input/geometry','eps');
DATA.inputs.SHEAR = h5readatt(filename,'/data/input/geometry','shear'); DATA.inputs.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
...@@ -38,10 +37,12 @@ function [DATA] = load_params(DATA,filename) ...@@ -38,10 +37,12 @@ function [DATA] = load_params(DATA,filename)
DATA.inputs.MUz = h5readatt(filename,'/data/input/model','mu_z'); DATA.inputs.MUz = h5readatt(filename,'/data/input/model','mu_z');
DATA.inputs.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY'); DATA.inputs.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
DATA.inputs.BETA = h5readatt(filename,'/data/input/model','beta'); DATA.inputs.BETA = h5readatt(filename,'/data/input/model','beta');
DATA.inputs.TAU_E = h5readatt(filename,'/data/input/model','tau_e'); DATA.inputs.TAU_I = h5readatt(filename,'/data/input/model','tau_i');
DATA.inputs.HYP_V = h5readatt(filename,'/data/input/model','HYP_V'); DATA.inputs.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
DATA.inputs.K_cB = h5readatt(filename,'/data/input/model','k_cB'); DATA.inputs.K_cB = h5readatt(filename,'/data/input/model','k_cB');
DATA.inputs.K_gB = h5readatt(filename,'/data/input/model','k_gB'); DATA.inputs.K_gB = h5readatt(filename,'/data/input/model','k_gB');
DATA.inputs.ADIAB_E = h5readatt(filename,'/data/input/model','ADIAB_E') == 'y';
DATA.inputs.ADIAB_I = h5readatt(filename,'/data/input/model','ADIAB_I') == 'y';
DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y'; DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y'; DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
...@@ -56,6 +57,13 @@ function [DATA] = load_params(DATA,filename) ...@@ -56,6 +57,13 @@ function [DATA] = load_params(DATA,filename)
DATA.inputs.K_N = zeros(1,DATA.inputs.Na); DATA.inputs.K_N = zeros(1,DATA.inputs.Na);
DATA.inputs.K_T = zeros(1,DATA.inputs.Na); DATA.inputs.K_T = zeros(1,DATA.inputs.Na);
spnames = {'ions','electrons'}; spnames = {'ions','electrons'};
if(DATA.inputs.ADIAB_E)
spnames = {spnames{1}};
end
if(DATA.inputs.ADIAB_I)
spnames = {spnames{2}};
end
for ia=1:DATA.inputs.Na for ia=1:DATA.inputs.Na
spdata = ['/data/input/',spnames{ia}]; spdata = ['/data/input/',spnames{ia}];
DATA.inputs.sigma(ia) = h5readatt(filename,spdata,'sigma'); DATA.inputs.sigma(ia) = h5readatt(filename,spdata,'sigma');
......
matlab/plot/plot_ballooning.m
+ 21
− 16
View file @ 79443f28
...@@ -2,20 +2,19 @@ function [FIG] = plot_ballooning(data,options) ...@@ -2,20 +2,19 @@ function [FIG] = plot_ballooning(data,options)
FIG.fig = figure; iplot = 1; FIG.fig = figure; iplot = 1;
[~,it0] = min(abs(data.Ts3D - options.time_2_plot(1))); [~,it0] = min(abs(data.Ts3D - options.time_2_plot(1)));
[~,it1] = min(abs(data.Ts3D - options.time_2_plot(end))); [~,it1] = min(abs(data.Ts3D - options.time_2_plot(end)));
[~,ikyarray] = min(abs(data.ky - options.kymodes)); [~,ikyarray] = min(abs(data.grids.ky - options.kymodes));
% phi_real=mean(real(data.PHI(:,:,:,it0:it1)),4); ikyarray = unique(ikyarray);
% phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
phi_real=real(data.PHI(:,:,:,it1)); phi_real=real(data.PHI(:,:,:,it1));
phi_imag=imag(data.PHI(:,:,:,it1)); phi_imag=imag(data.PHI(:,:,:,it1));
ncol = 1; ncol = 1;
if data.BETA > 0 if data.inputs.BETA > 0
psi_real=real(data.PSI(:,:,:,it1)); psi_real=real(data.PSI(:,:,:,it1));
psi_imag=imag(data.PSI(:,:,:,it1)); psi_imag=imag(data.PSI(:,:,:,it1));
ncol = 2; ncol = 2;
end end
% Apply ballooning transform % Apply ballooning transform
if(data.Nkx > 1) if(data.grids.Nkx > 1)
nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2)); nexc = round(data.grids.ky(2)*data.inputs.SHEAR*2*pi/data.grids.kx(2));
else else
nexc = 1; nexc = 1;
end end
...@@ -34,7 +33,7 @@ function [FIG] = plot_ballooning(data,options) ...@@ -34,7 +33,7 @@ function [FIG] = plot_ballooning(data,options)
ordered_ikx = [(Np_+1):Nkx 1:Np_]; ordered_ikx = [(Np_+1):Nkx 1:Np_];
end end
try try
idx=data.kx./data.kx(2); idx=data.grids.kx./data.grids.kx(2);
catch catch
idx=0; idx=0;
end end
...@@ -54,11 +53,11 @@ function [FIG] = plot_ballooning(data,options) ...@@ -54,11 +53,11 @@ function [FIG] = plot_ballooning(data,options)
end end
% Define ballooning angle % Define ballooning angle
coordz = data.z; coordz = data.grids.z;
for i_ =1: Nkx for i_ =1: Nkx
for iz=1:dims(3) for iz=1:dims(3)
ii = dims(3)*(i_-1) + iz; ii = dims(3)*(i_-1) + iz;
b_angle(ii) =coordz(iz) + 2*pi*idx(i_)./is; b_angle(ii) =coordz(iz) + 2*pi*data.grids.Npol*idx(i_)./is;
end end
end end
...@@ -73,6 +72,7 @@ function [FIG] = plot_ballooning(data,options) ...@@ -73,6 +72,7 @@ function [FIG] = plot_ballooning(data,options)
phib_norm = phib / normalization; phib_norm = phib / normalization;
phib_real_norm = real( phib_norm); phib_real_norm = real( phib_norm);
phib_imag_norm = imag( phib_norm); phib_imag_norm = imag( phib_norm);
%
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+1) subplot(numel(ikyarray),ncol,ncol*(iplot-1)+1)
plot(b_angle/pi, phib_real_norm,'-b'); hold on; plot(b_angle/pi, phib_real_norm,'-b'); hold on;
plot(b_angle/pi, phib_imag_norm ,'-r'); plot(b_angle/pi, phib_imag_norm ,'-r');
...@@ -80,18 +80,20 @@ function [FIG] = plot_ballooning(data,options) ...@@ -80,18 +80,20 @@ function [FIG] = plot_ballooning(data,options)
legend('real','imag','norm') legend('real','imag','norm')
xlabel('$\chi / \pi$') xlabel('$\chi / \pi$')
ylabel('$\phi_B (\chi)$'); ylabel('$\phi_B (\chi)$');
title(['$k_y=',sprintf('%2.2f',data.ky(iky)),... title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']); ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
for i_ = 2*[1:data.grids.Nkx]-(data.grids.Nkx+1)
if data.BETA > 0 xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
if data.inputs.BETA > 0
psib_real = zeros( Nkx*dims(3) ,1); psib_real = zeros( Nkx*dims(3) ,1);
psib_imag = psib_real; psib_imag = psib_real;
for i_ =1:Nkx for i_ =1:Nkx
start_ = (i_ -1)*dims(3) +1; start_ = (i_-1)*dims(3) +1;
end_ = i_*dims(3); end_ = i_*dims(3);
ikx = ordered_ikx(i_); ikx = ordered_ikx(i_);
psib_real(start_:end_) = psi_real(iky,ikx,:); psib_real(start_:end_) = psi_real(iky,ikx,:);
psib_imag(start_:end_) = psi_imag(iky,ikx,:); psib_imag(start_:end_) = psi_imag(iky,ikx,:);
end end
psib = psib_real(:) + 1i * psib_imag(:); psib = psib_real(:) + 1i * psib_imag(:);
psib_norm = psib / normalization; psib_norm = psib / normalization;
...@@ -102,9 +104,12 @@ function [FIG] = plot_ballooning(data,options) ...@@ -102,9 +104,12 @@ function [FIG] = plot_ballooning(data,options)
plot(b_angle/pi, psib_imag_norm ,'-r'); plot(b_angle/pi, psib_imag_norm ,'-r');
plot(b_angle/pi, abs(psib_norm),'-k'); plot(b_angle/pi, abs(psib_norm),'-k');
legend('real','imag','norm') legend('real','imag','norm')
for i_ = 2*[1:data.grids.Nkx]-(data.grids.Nkx+1)
xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
xlabel('$\chi / \pi$') xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$'); ylabel('$\psi_B (\chi)$');
title(['$k_y=',sprintf('%2.2f',data.ky(iky)),... title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']); ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
end end
......
matlab/plot/plot_cosol_mat.m
+ 10
− 17
View file @ 79443f28
...@@ -108,9 +108,10 @@ mfns = {... ...@@ -108,9 +108,10 @@ mfns = {...
'/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';... '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';...
'/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';... '/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';...
'/home/ahoffman/gyacomo/iCa/gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';... '/home/ahoffman/gyacomo/iCa/gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';...
'/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';... '/home/ahoffman/gyacomo/iCa/gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';...
'/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';... % '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';...
'/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';... % '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
% '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5';... % '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5';...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5';... % '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5';...
}; };
...@@ -125,21 +126,13 @@ CONAME_A = {... ...@@ -125,21 +126,13 @@ CONAME_A = {...
'FC 4 2 NFLR 12'; ... 'FC 4 2 NFLR 12'; ...
'LD 10 5 NFLR 12'; ... 'LD 10 5 NFLR 12'; ...
'LD 11 7 NFLR 16'; ... 'LD 11 7 NFLR 16'; ...
'SG 11 7 NFLR 12, tau 1e-3'; ... 'LDii 16 9 NFLR 8'; ...
'SG 4 2 NFLR 5, tau 1e-3'; ... % 'SG 11 7 NFLR 12, tau 1e-3'; ...
'SG 4 2 NFLR 5, tau 1e-2'; ... % 'SG 4 2 NFLR 5, tau 1e-3'; ...
% 'SG 4 2 NFLR 5, tau 1e-2'; ...
% 'Hacked SG A';... % 'Hacked SG A';...
% 'Hacked SG B';... % 'Hacked SG B';...
}; };
TAU_A = [1;...
1;...
1;...
1;...
1;...
1e-3;...
1e-3;...
1e-2;...
];
figure figure
for j_ = 1:numel(mfns) for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name); mat_file_name = mfns{j_}; disp(mat_file_name);
...@@ -154,9 +147,9 @@ for j_ = 1:numel(mfns) ...@@ -154,9 +147,9 @@ for j_ = 1:numel(mfns)
wmax(idx_+1) = max(abs(imag(eig(MAT)))); wmax(idx_+1) = max(abs(imag(eig(MAT))));
end end
subplot(121) subplot(121)
plot(kp_a,gmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on; plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
subplot(122) subplot(122)
plot(kp_a,wmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on; plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
end end
subplot(121) subplot(121)
legend('show'); grid on; legend('show'); grid on;
......
matlab/plot/plot_metric.m
+ 15
− 14
View file @ 79443f28
...@@ -5,42 +5,43 @@ function [ fig, arrays ] = plot_metric( data, options ) ...@@ -5,42 +5,43 @@ function [ fig, arrays ] = plot_metric( data, options )
names = {'gxx','gxy','gxz','gyy','gyz','gzz',... names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
'hatB', 'dBdx', 'dBdy', 'dBdz',... 'hatB', 'dBdx', 'dBdy', 'dBdz',...
'Jacobian','hatR','hatZ','gradz_coeff'}; 'Jacobian','hatR','hatZ','gradz_coeff'};
geo_arrays = zeros(2,data.Nz,numel(names)); geo_arrays = zeros(data.grids.Nz,numel(names));
for i_ = 1:numel(names) for i_ = 1:numel(names)
namae = names{i_}; namae = names{i_};
geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])'; geo_arrays(:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end end
NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS; NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
if NPLOT > 0 if NPLOT > 0
fig = figure; fig = figure;
if options.SHOW_METRICS if options.SHOW_METRICS
subplot(311) subplot(3,NPLOT,1*NPLOT)
for i = 1:6 for i = 1:6
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry'); xlim([min(data.grids.z),max(data.grids.z)]); legend('show'); title('GYACOMO geometry');
subplot(312) subplot(3,NPLOT,2*NPLOT)
for i = 7:10 for i = 7:10
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(data.z),max(data.z)]); legend('show'); xlim([min(data.grids.z),max(data.grids.z)]); legend('show');
subplot(313) subplot(3,NPLOT,3*NPLOT)
for i = 11:14 for i = 11:14
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(data.z),max(data.z)]); legend('show'); xlim([min(data.grids.z),max(data.grids.z)]); legend('show');
end end
if options.SHOW_FLUXSURF if options.SHOW_FLUXSURF
R = [geo_arrays(1,:,12),geo_arrays(1,1,12)]; subplot(1,2,1)
Z = [geo_arrays(1,:,13),geo_arrays(1,1,13)]; R = [geo_arrays(:,12);geo_arrays(1,12)];
Z = [geo_arrays(:,13);geo_arrays(1,13)];
plot(R,Z); plot(R,Z);
axis equal axis equal
end end
end end
%outputs %outputs
arrays = squeeze(geo_arrays(1,:,:)); arrays = squeeze(geo_arrays(:,:));
end end
matlab/plot/show_geometry.m
+ 12
− 12
View file @ 79443f28
function [ FIGURE ] = show_geometry(DATA,OPTIONS) function [ FIGURE ] = show_geometry(DATA,OPTIONS)
% filtering Z pinch and torus % filtering Z pinch and torus
if DATA.Nz > 1 % Torus flux tube geometry if DATA.grids.Nz > 1 % Torus flux tube geometry
Nturns = floor(abs(DATA.z(1)/pi)); Nturns = floor(abs(DATA.grids.z(1)/pi));
Nptor = ceil(DATA.Nz*2/Nturns); Tgeom = 1; Nptor = ceil(DATA.grids.Nz*2/Nturns); Tgeom = 1;
a = DATA.EPS; % Torus minor radius a = DATA.EPS; % Torus minor radius
R = 1.; % Torus major radius R = 1.; % Torus major radius
q = (DATA.Nz>1)*DATA.Q0; % Flux tube safety factor q = (DATA.grids.Nz>1)*DATA.inputs.Q0; % Flux tube safety factor
theta = linspace(-pi, pi, Nptor) ; % Poloidal angle theta = linspace(-pi, pi, Nptor) ; % Poloidal angle
phi = linspace(0., 2.*pi, Nptor) ; % Toroidal angle phi = linspace(0., 2.*pi, Nptor) ; % Toroidal angle
[t, p] = meshgrid(phi, theta); [t, p] = meshgrid(phi, theta);
x_tor = (R + a.*cos(p)) .* cos(t); x_tor = (R + a.*cos(p)) .* cos(t);
y_tor = (R + a.*cos(p)) .* sin(t); y_tor = (R + a.*cos(p)) .* sin(t);
z_tor = a.*sin(p); z_tor = a.*sin(p);
DIMENSIONS = [50 50 1200 600]; DIMENSIONS = [600 600 1200 600];
else % Zpinch geometry else % Zpinch geometry
Nturns = 0.1; Tgeom = 0; Nturns = 0.1; Tgeom = 0;
a = 0.7; % Torus minor radius a = 0.7; % Torus minor radius
...@@ -48,7 +48,7 @@ yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z) ...@@ -48,7 +48,7 @@ yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z)
yvec= @(z) [yX(z); yY(z); yZ(z)]; yvec= @(z) [yX(z); yY(z); yZ(z)];
% Plot high res field line % Plot high res field line
% Planes plot % Planes plot
theta = DATA.z ; % Poloidal angle theta = DATA.grids.z ; % Poloidal angle
% Plot x,y,bvec at these points % Plot x,y,bvec at these points
scale = 0.10; scale = 0.10;
...@@ -56,11 +56,11 @@ scale = 0.10; ...@@ -56,11 +56,11 @@ scale = 0.10;
OPTIONS.POLARPLOT = 0; OPTIONS.POLARPLOT = 0;
OPTIONS.PLAN = 'xy'; OPTIONS.PLAN = 'xy';
r_o_R = DATA.rho_o_R; r_o_R = DATA.rho_o_R;
[X,Y] = meshgrid(r_o_R*DATA.x,r_o_R*DATA.y); [X,Y] = meshgrid(r_o_R*DATA.grids.x,r_o_R*DATA.grids.y);
max_ = 0; max_ = 0;
FIELDS = zeros(DATA.Ny,DATA.Nx,DATA.Nz); FIELDS = zeros(DATA.grids.Ny,DATA.grids.Nx,DATA.grids.Nz);
FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.PARAMS]; set(gcf, 'Position', DIMENSIONS) FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.params_string]; set(gcf, 'Position', DIMENSIONS)
for it_ = 1:numel(OPTIONS.TIME) for it_ = 1:numel(OPTIONS.TIME)
subplot(1,numel(OPTIONS.TIME),it_) subplot(1,numel(OPTIONS.TIME),it_)
%plot magnetic geometry %plot magnetic geometry
...@@ -86,7 +86,7 @@ subplot(1,numel(OPTIONS.TIME),it_) ...@@ -86,7 +86,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
end end
end end
%plot vector basis %plot vector basis
theta = DATA.z ; % Poloidal angle theta = DATA.grids.z ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on; plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r'); quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r');
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g'); quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
...@@ -101,7 +101,7 @@ subplot(1,numel(OPTIONS.TIME),it_) ...@@ -101,7 +101,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
if (tmp > max_) max_ = tmp; if (tmp > max_) max_ = tmp;
end end
for iz = OPTIONS.PLANES for iz = OPTIONS.PLANES
z_ = DATA.z(iz); z_ = DATA.grids.z(iz);
Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_); Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_); Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_); Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
...@@ -110,7 +110,7 @@ subplot(1,numel(OPTIONS.TIME),it_) ...@@ -110,7 +110,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
colormap(bluewhitered); colormap(bluewhitered);
% caxis([-1,1]); % caxis([-1,1]);
end end
if DATA.Nz == 1 if DATA.grids.Nz == 1
xlim([0.65 0.75]); xlim([0.65 0.75]);
ylim([-0.1 0.1]); ylim([-0.1 0.1]);
zlim([-0.2 0.2]); zlim([-0.2 0.2]);
......
matlab/setup.m
+ 16
− 5
View file @ 79443f28
...@@ -15,17 +15,28 @@ GEOM.geom = ['''',GEOMETRY,'''']; ...@@ -15,17 +15,28 @@ GEOM.geom = ['''',GEOMETRY,''''];
GEOM.q0 = Q0; % q factor GEOM.q0 = Q0; % q factor
GEOM.shear = SHEAR; % shear GEOM.shear = SHEAR; % shear
GEOM.eps = EPS; % inverse aspect ratio GEOM.eps = EPS; % inverse aspect ratio
GEOM.kappa = KAPPA; % elongation GEOM.kappa = KAPPA; % elongation
GEOM.delta = DELTA; % triangularity GEOM.s_kappa = S_KAPPA;
GEOM.zeta = ZETA; % squareness GEOM.delta = DELTA; % triangularity
GEOM.s_delta = S_DELTA;
GEOM.zeta = ZETA; % squareness
GEOM.s_zeta = S_ZETA;
GEOM.parallel_bc = ['''',PARALLEL_BC,'''']; GEOM.parallel_bc = ['''',PARALLEL_BC,''''];
GEOM.shift_y = SHIFT_Y; GEOM.shift_y = SHIFT_Y;
GEOM.Npol = NPOL; GEOM.Npol = NPOL;
if PB_PHASE; GEOM.PB_PHASE = '.true.'; else; GEOM.PB_PHASE = '.false.';end;
% Model parameters % Model parameters
MODEL.LINEARITY = ['''',LINEARITY,'''']; MODEL.LINEARITY = ['''',LINEARITY,''''];
try
RM_LD_T_EQ;
catch
RM_LD_T_EQ = 0;
end
if RM_LD_T_EQ; MODEL.RM_LD_T_EQ = '.true.'; else; MODEL.RM_LD_T_EQ = '.false.'; end; if RM_LD_T_EQ; MODEL.RM_LD_T_EQ = '.true.'; else; MODEL.RM_LD_T_EQ = '.false.'; end;
MODEL.Na = NA; MODEL.Na = NA;
if ADIAB_E; MODEL.ADIAB_E = '.true.'; else; MODEL.ADIAB_E = '.false.';end; if ADIAB_E; MODEL.ADIAB_E = '.true.'; else; MODEL.ADIAB_E = '.false.';end;
if ADIAB_I; MODEL.ADIAB_I = '.true.'; else; MODEL.ADIAB_I = '.false.';end;
if MHD_PD; MODEL.MHD_PD = '.true.'; else; MODEL.MHD_PD = '.false.';end;
MODEL.beta = BETA; MODEL.beta = BETA;
MODEL.mu_x = MU_X; MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y; MODEL.mu_y = MU_Y;
...@@ -44,7 +55,6 @@ MODEL.sigma_i = 1.0; ...@@ -44,7 +55,6 @@ MODEL.sigma_i = 1.0;
% charge q_a/e % charge q_a/e
MODEL.q_e =-1.0; MODEL.q_e =-1.0;
MODEL.q_i = 1.0; MODEL.q_i = 1.0;
if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x % gradients L_perp/L_x
MODEL.K_Ni = K_Ni; MODEL.K_Ni = K_Ni;
MODEL.K_Ne = K_Ne; MODEL.K_Ne = K_Ne;
...@@ -71,7 +81,8 @@ switch CO ...@@ -71,7 +81,8 @@ switch CO
case 'LR' case 'LR'
COLL.mat_file = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; COLL.mat_file = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
case 'LD' case 'LD'
COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5'; % COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
COLL.mat_file = 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';
% COLL.mat_file = 'gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5'; % COLL.mat_file = 'gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';
% COLL.mat_file = 'gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'; % COLL.mat_file = 'gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
% COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5'; % COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';
......
matlab/write_fort90.m
+ 20
− 12
View file @ 79443f28
...@@ -29,11 +29,15 @@ fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']); ...@@ -29,11 +29,15 @@ fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']);
fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']); fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']);
fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']); fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']); fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
fprintf(fid,['s_kappa = ', num2str(GEOM.s_kappa),'\n']);
fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']); fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
fprintf(fid,['s_delta = ', num2str(GEOM.s_delta),'\n']);
fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']); fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
fprintf(fid,['s_zeta = ', num2str(GEOM.s_zeta),'\n']);
fprintf(fid,[' parallel_bc = ', GEOM.parallel_bc,'\n']); fprintf(fid,[' parallel_bc = ', GEOM.parallel_bc,'\n']);
fprintf(fid,[' shift_y = ', num2str(GEOM.shift_y),'\n']); fprintf(fid,[' shift_y = ', num2str(GEOM.shift_y),'\n']);
fprintf(fid,[' Npol = ', num2str(GEOM.Npol),'\n']); fprintf(fid,[' Npol = ', num2str(GEOM.Npol),'\n']);
fprintf(fid,[' PB_PHASE= ', GEOM.PB_PHASE,'\n']);
fprintf(fid,'/\n'); fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n'); fprintf(fid,'&OUTPUT_PAR\n');
...@@ -69,6 +73,9 @@ fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']); ...@@ -69,6 +73,9 @@ fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']); fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']); fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
fprintf(fid,[' ADIAB_E = ', MODEL.ADIAB_E,'\n']); fprintf(fid,[' ADIAB_E = ', MODEL.ADIAB_E,'\n']);
fprintf(fid,[' ADIAB_I = ', MODEL.ADIAB_I,'\n']);
fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
fprintf(fid,[' MHD_PD = ', MODEL.MHD_PD,'\n']);
fprintf(fid,'/\n'); fprintf(fid,'/\n');
fprintf(fid,'&CLOSURE_PAR\n'); fprintf(fid,'&CLOSURE_PAR\n');
...@@ -78,18 +85,19 @@ fprintf(fid,[' nonlinear_closure=',CLOSURE.nonlinear_closure,'\n']); ...@@ -78,18 +85,19 @@ fprintf(fid,[' nonlinear_closure=',CLOSURE.nonlinear_closure,'\n']);
fprintf(fid,[' nmax =',num2str(CLOSURE.nmax),'\n']); fprintf(fid,[' nmax =',num2str(CLOSURE.nmax),'\n']);
fprintf(fid,'/\n'); fprintf(fid,'/\n');
fprintf(fid,'&SPECIES\n'); if(strcmp(MODEL.ADIAB_I,'.false.'))
fprintf(fid, ' name_ = ions \n'); fprintf(fid,'&SPECIES\n');
fprintf(fid,[' tau_ = ', num2str(MODEL.tau_i),'\n']); fprintf(fid, ' name_ = ions \n');
fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_i),'\n']); fprintf(fid,[' tau_ = ', num2str(MODEL.tau_i),'\n']);
fprintf(fid,[' q_ = ', num2str(MODEL.q_i),'\n']); fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ni),'\n']); fprintf(fid,[' q_ = ', num2str(MODEL.q_i),'\n']);
fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Ti),'\n']); fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ni),'\n']);
fprintf(fid,'/\n'); fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Ti),'\n']);
fprintf(fid,'/\n');
if(MODEL.Na > 1) end
fprintf(fid,'&SPECIES\n'); if(strcmp(MODEL.ADIAB_E,'.false.'))
fprintf(fid, ' name_ = electrons'); fprintf(fid,'&SPECIES\n');
fprintf(fid, ' name_ = electrons \n');
fprintf(fid,[' tau_ = ', num2str(MODEL.tau_e),'\n']); fprintf(fid,[' tau_ = ', num2str(MODEL.tau_e),'\n']);
fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_e),'\n']); fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' q_ = ', num2str(MODEL.q_e),'\n']); fprintf(fid,[' q_ = ', num2str(MODEL.q_e),'\n']);
......
wk/load_metadata_scan.m → wk/CBC_load_metadata_scan.m
+ 24
− 12
View file @ 79443f28
...@@ -60,15 +60,26 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0' ...@@ -60,15 +60,26 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% 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.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'; % 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_64_DGDK_0.001_kT_5.3.mat';
%% ky pj scan %% 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_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_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_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_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'; % 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_P_4_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_30_DGDK_0.001.mat'; % datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_4_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_4_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_10_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
%% Chose if we filter gamma>0.05 %% Chose if we filter gamma>0.05
FILTERGAMMA = 1; FILTERGAMMA = 1;
...@@ -104,14 +115,15 @@ if 1 ...@@ -104,14 +115,15 @@ if 1
figure figure
colors_ = jet(numel(d.s2)); colors_ = jet(numel(d.s2));
for i = 1:numel(d.s2) for i = 1:numel(d.s2)
% plot(d.s1,d.data(:,i),'s-',... % plot(d.s1,d.data(:,i),'s-',...
% 'LineWidth',2.0,... plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],... 'LineWidth',2.0,...
% 'color',colors_(i,:)); 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',... 'color',colors_(i,:));
'LineWidth',2.0,... % errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],... % 'LineWidth',2.0,...
'color',colors_(i,:)); % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
hold on; hold on;
end end
xlabel(d.s1name); ylabel(d.dname);title(d.title); xlabel(d.s1name); ylabel(d.dname);title(d.title);
......
wk/ETG_adiab_i.m 0 → 100644
+ 162
− 0
View file @ 79443f28
%% QUICK RUN SCRIPT
% This script creates a directory in /results and runs a simulation directly
% from the Matlab framework. It is meant to run only small problems in linear
% for benchmarking and debugging purposes since it makes Matlab "busy".
%% Set up the paths for the necessary Matlab modules
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters
SIMID = 'ETG_adiab_i'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 0.0; % Collision frequency
TAU = 1; % e/i temperature ratio
K_Ne = 0; % ele Density '''
K_Te = 10; % ele Temperature '''
K_Ni = 0; % ion Density gradient drive
K_Ti = 0; % ion Temperature '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 1; % number of kinetic species
ADIAB_E = 0; % adiabatic electron model
ADIAB_I = 1; % adiabatic ion model
BETA = 1e-5; % electron plasma beta
MHD_PD = 0; % MHD pressure drift
%% Set up grid parameters
P = 2;
J = P/2;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 2; % real space x-gridpoints
NY = 16; % real space y-gridpoints
LX = 2*pi/0.3; % Size of the squared frequency domain in x direction
LY = 2*pi/4.0; % Size of the squared frequency domain in y direction
NZ = 16; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
% GEOMETRY= 'miller';
GEOMETRY= 'z-pinch';
EPS = 1.0; % inverse aspect ratio
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
KAPPA = 1.0; % elongation
S_KAPPA = 0.0;
DELTA = 0.0; % triangularity
S_DELTA = 0.0;
ZETA = 0.0; % squareness
S_ZETA = 0.0;
PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
SHIFT_Y = 0.0; % Shift in the periodic BC in z
NPOL = 1; % Number of poloidal turns
PB_PHASE= 0;
%% TIME PARAMETERS
TMAX = 10; % Maximal time unit
DT = 1e-3; % Time step
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 3; % 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)
CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
KERN = 0; % Kernel model (0 : GK)
INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
%% OUTPUTS
W_DOUBLE = 1; % Output flag for double moments
W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
W_HF = 1; % Output flag for high-frequency potential energy
W_PHI = 1; % Output flag for potential
W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
W_DENS = 1; % Output flag for total density
W_TEMP = 1; % Output flag for temperature
W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% UNUSED PARAMETERS
% These parameters are usually not to play with in linear runs
MU = 0.1; % Hyperdiffusivity coefficient
MU_X = MU; % Hyperdiffusivity coefficient in x direction
MU_Y = MU; % Hyperdiffusivity coefficient in y direction
N_HD = 4; % Degree of spatial-hyperdiffusivity
MU_Z = 10.0; % Hyperdiffusivity coefficient in z direction
HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
LAMBDAD = 0.0; % Lambda Debye
NOISE0 = 0.0e-5; % Initial noise amplitude
BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
% RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
% Load outputs from jobnummin up to jobnummax
J0 = 0; J1 = 0;
data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
if 1 % Activate or not
%% plot mode evolution and growth rates
% Load phi
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
options.NORMALIZED = 0;
options.TIME = data.Ts3D;
% Time window to measure the growth of kx/ky modes
options.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*data.Ts3D(end);
options.NMA = 1; % Set NMA option to 1
options.NMODES = 999; % Set how much modes we study
options.iz = 'avg'; % Compressing z
options.ik = 1; %
options.fftz.flag = 0; % Set fftz.flag option to 0
fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
end
wk/KBM_ky_PJ_scan.m 0 → 100644
+ 218
− 0
View file @ 79443f28
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 = 'lin_KBM'; % Name of the simulation
RERUN = 1; % rerun if the data does not exist
RUN = 1;
K_Ne = 3; % ele Density '''
K_Te = 4.5; % ele Temperature '''
K_Ni = 3; % ion Density gradient drive
K_Ti = 8; % ion Temperature '''
P_a = [2 4 8 16];
% P_a = 10;
ky_a= 0.05:0.1:0.75;
% ky_a = 0.05;
% collision setting
CO = 'DG';
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
NU = 1e-2;
% model
NA = 2;
BETA = 0.03; % electron plasma beta
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 5e-3;
TMAX = 15;
% 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)
%% GRID PARAMETERS
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = P/2; % Laguerre "
NX = 12; % 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 = 0.5; % 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);
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)>5
% 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),'_be_',num2str(BETA),'.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/KBM_load_metadata_scan.m 0 → 100644
+ 171
− 0
View file @ 79443f28
% Metadata path
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';
datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 0;
%% Load data
fname = ['../results/',datafname];
d = load(fname);
if FILTERGAMMA
d.data = d.data.*(d.data>0.025);
d.err = d.err.*(d.data>0.025);
end
if 0
%% Pcolor of the peak
figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
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);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
colormap(jet)
colormap(bluewhitered)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 1
%% Scan along first dimension
figure
colors_ = jet(numel(d.s2));
for i = 1:numel(d.s2)
% plot(d.s1,d.data(:,i),'s-',...
plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
% errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
colormap(colors_);
clb = colorbar;
% 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 =1.5:numel(d.s2)+1.5;
clb.TickLabels=d.s2;
clb.Label.String = d.s2name;
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 0
%% Scan along second dimension
figure
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
plot(d.s2,d.data(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
% 'color',colors_(i,:));
hold on;
end
xlabel(d.s2name); ylabel(d.dname);title(d.title);
xlim([d.s2(1) d.s2(end)]);
colormap(jet(numel(d.s1)));
clb = colorbar;
caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 0
%% Convergence analysis
figure
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
% target_ = d.data(i,end);
target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
if target_ > 0
eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
semilogy(d.s2(1:end-1),eps_,'-s',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
hold on;
end
end
xlabel(d.s2name); ylabel('$\epsilon_r$');title(d.title);
xlim([d.s2(1) d.s2(end)]);
colormap(colors_);
clb = colorbar;
caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
grid on;
end
if 0
%% Pcolor of the error
figure; i_ = 0;
% 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.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_));
% eps_ = log(abs(target_ - d.data)/abs(target_));
eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
sign_ = 1;%sign(d.data - target_);
eps_ = d.data;
for i = 1:numel(d.s1)
for j = 1:numel(d.s2)
% target_ = d.data(i,end);
% target_ = d.data(i,end);
% target_ = d.data(end,j);
eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
if target_ > 0
% 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
[XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
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);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
% colormap(jet)
colormap(bluewhitered)
% caxis([-10, 0]);
clb=colorbar;
clb.Label.String = '$\log(\epsilon_r)$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
\ No newline at end of file
wk/analysis_gene.m
+ 8
− 8
View file @ 79443f28
...@@ -60,9 +60,9 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add ...@@ -60,9 +60,9 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% 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/30x16x128x64x24/kT_4.0/';
folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/'; % 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/Miller_KT_6.96_128x64x24x16x8/output/';
% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/'; % folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
% debug ? shearless % debug ? shearless
...@@ -176,21 +176,21 @@ names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',... ...@@ -176,21 +176,21 @@ names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
figure; figure;
subplot(311) subplot(311)
for i = 1:6 for i = 1:6
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry'); xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show'); title('GENE geometry');
subplot(312) subplot(312)
for i = 7:10 for i = 7:10
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show');
subplot(313) subplot(313)
for i = [11 12 14] for i = [11 12 14]
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show');
end end
......
wk/analysis_gyacomo.m
+ 2
− 2
View file @ 79443f28
...@@ -137,7 +137,7 @@ options.PLT_FTUBE = 0; ...@@ -137,7 +137,7 @@ options.PLT_FTUBE = 0;
data.EPS = 0.4; data.EPS = 0.4;
data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options); fig = show_geometry(data,options);
save_figure(data,fig,'.png') % save_figure(data,fig,'.png')
end end
if 0 if 0
...@@ -231,7 +231,7 @@ end ...@@ -231,7 +231,7 @@ end
if 0 if 0
%% Metric infos %% Metric infos
options.SHOW_FLUXSURF = 0; options.SHOW_FLUXSURF = 1;
options.SHOW_METRICS = 1; options.SHOW_METRICS = 1;
[fig, geo_arrays] = plot_metric(data,options); [fig, geo_arrays] = plot_metric(data,options);
end end
......
wk/benchmark_and_scan_scripts/Lin_1999_fig3_heat_diff_vs_time.txt 0 → 100644
+ 85
− 0
View file @ 79443f28
22.939632545931715, 0.0028421861631937606
46.876640419947535, 0.0015479763924537426
73.80577427821521, 0.0016458578036862015
91.75853018372703, 0.016904482243587093
106.71916010498688, 0.043201954728034675
117.69028871391077, 0.07639100361073647
123.67454068241472, 0.11646800365423937
129.65879265091866, 0.16897594292426155
131.65354330708664, 0.218706950305245
139.63254593175856, 0.32094589695625053
143.62204724409452, 0.4217891271878308
151.60104986876644, 0.5585584605791679
155.5905511811024, 0.6994569394295327
162.57217847769027, 0.8541784486883891
169.5538057742782, 0.9453640463450356
177.53280839895018, 0.7313046504546048
182.51968503937013, 0.5172343788518148
185.5118110236221, 0.4426596192050579
196.482939632546, 0.365351430518699
209.44881889763792, 0.2797631994895665
225.4068241469817, 0.23700352373080436
242.36220472440954, 0.21358448978408084
256.3254593175854, 0.1597678397935065
264.3044619422573, 0.12250402401357285
285.249343832021, 0.09357462913820858
308.1889763779529, 0.06050883832891052
319.1601049868768, 0.06331114688012063
330.13123359580055, 0.07301953277939699
346.0892388451445, 0.06617145923057977
357.0603674540684, 0.05239918214643069
376.0104986876642, 0.05661170806687843
390.97112860892406, 0.05804730209828746
415.9055118110238, 0.05123185568654742
450.81364829396335, 0.04721509258856482
480.73490813648306, 0.05699235799944913
503.67454068241494, 0.07917518597468143
531.6010498687665, 0.08756398544104638
553.5433070866145, 0.08902495613462691
580.4724409448822, 0.11674714693812438
600.4199475065618, 0.13477545279215797
621.3648293963256, 0.14452009106596486
636.3254593175855, 0.13214353040124127
651.2860892388453, 0.16120343382491553
661.2598425196852, 0.18748278012209796
676.2204724409451, 0.21654268354577222
687.191601049869, 0.27597482635112613
690.1837270341209, 0.31189730427343
694.1732283464569, 0.3464421919635735
700.1574803149608, 0.3727070373109439
702.1522309711288, 0.3768579341946898
717.1128608923887, 0.3755310972868723
723.0971128608926, 0.3548346166673918
727.0866141732286, 0.3216999463464858
736.0629921259845, 0.27062760110787254
745.0393700787404, 0.24441713432229806
758.0052493438322, 0.18507199721581769
768.9763779527561, 0.1409129798001768
777.952755905512, 0.10779643566653618
785.931758530184, 0.08020112817389546
792.9133858267719, 0.0636519192007079
799.8950131233598, 0.0526275721059728
820.8398950131236, 0.04856005568364741
844.7769028871394, 0.04864706138252073
855.7480314960633, 0.044543292585664584
872.7034120734911, 0.047367352561592746
897.6377952755906, 0.06403256913327837
914.5931758530187, 0.07376270645727301
929.5538057742783, 0.07934194689752161
955.4855643044623, 0.08219863401052785
969.4488188976381, 0.08639303374371021
983.412073490814, 0.09611229535534571
989.396325459318, 0.09751526224967744
1008.3464566929135, 0.1031090036397384
1027.2965879265096, 0.11422760690825251
1037.2703412073495, 0.1253135830396891
1047.2440944881894, 0.13087469729267276
1061.2073490813652, 0.14888125172198774
1068.1889763779532, 0.16686242948913155
1076.167979002625, 0.1820848015545018
1084.1469816272968, 0.19868838908948538
1100.1049868766409, 0.22913313322022588
1113.0708661417327, 0.24713606241208796
1116.0629921259845, 0.24852815359406044
1126.0367454068244, 0.2416583286205246
1134.0157480314965, 0.2237315294151767
wk/benchmark_and_scan_scripts/old/fort_00.90
+ 31
− 23
View file @ 79443f28
&BASIC &BASIC
nrun = 100000000 nrun = 100000000
dt = 0.02 dt = 0.005
tmax = 30 tmax = 15
maxruntime = 356400 maxruntime = 356400
job2load = -1 job2load = -1
/ /
&GRID &GRID
pmax = 60 pmax = 4
jmax = 30 jmax = 2
Nx = 8 Nx = 9
Lx = 7.854 Lx = 62.8319
Ny = 2 Ny = 2
Ly = 20.944 Ly = 25.1327
Nz = 24 Nz = 16
SG = .false. SG = .false.
Nexc = 1 Nexc = 1
/ /
&GEOMETRY &GEOMETRY
geom = 's-alpha' geom = 'miller'
q0 = 1.4 q0 = 1.4
shear = 0.8 shear = 0.8
eps = 0.18 eps = 0.18
...@@ -34,31 +34,32 @@ ...@@ -34,31 +34,32 @@
dtsave_2d = -1 dtsave_2d = -1
dtsave_3d = 2 dtsave_3d = 2
dtsave_5d = 100 dtsave_5d = 100
write_doubleprecision = .false. write_doubleprecision = .true.
write_gamma = .true. write_gamma = .true.
write_hf = .true. write_hf = .true.
write_phi = .true. write_phi = .true.
write_Na00 = .true. write_Na00 = .true.
write_Napj = .false. write_Napj = .true.
write_dens = .false. write_dens = .true.
write_temp = .true. write_temp = .true.
/ /
&MODEL_PAR &MODEL_PAR
LINEARITY = 'linear' LINEARITY = 'linear'
Na = 1 RM_LD_T_EQ= .false.
Na = 2
mu_x = 0 mu_x = 0
mu_y = 0 mu_y = 0
N_HD = 4 N_HD = 4
mu_z = 1 mu_z = 2
HYP_V = 'hypcoll' HYP_V = 'hypcoll'
mu_p = 0 mu_p = 0
mu_j = 0 mu_j = 0
nu = 0.1 nu = 0.005
k_gB = 1 k_gB = 1
k_cB = 1 k_cB = 1
lambdaD = 0 lambdaD = 0
beta = 0.0001 beta = 0.03
ADIAB_E = .true. ADIAB_E = .false.
/ /
&CLOSURE_PAR &CLOSURE_PAR
hierarchy_closure='truncation' hierarchy_closure='truncation'
...@@ -71,20 +72,27 @@ ...@@ -71,20 +72,27 @@
tau_ = 1 tau_ = 1
sigma_ = 1 sigma_ = 1
q_ = 1 q_ = 1
K_N_ = 2.22 K_N_ = 3
K_T_ = 3.5 K_T_ = 8
/
&SPECIES
name_ = electrons tau_ = 1
sigma_ = 0.051962
q_ = -1
K_N_ = 3
K_T_ = 4.5
/ /
&COLLISION_PAR &COLLISION_PAR
collision_model = 'DG' collision_model = 'DG'
GK_CO = .false. GK_CO = .false.
INTERSPECIES = .true. INTERSPECIES = .true.
mat_file = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/oiCa/null' mat_file = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/oiCa/null'
collision_kcut = 1.75 collision_kcut = 1
/ /
&INITIAL_CON &INITIAL_CON
INIT_OPT = 'phi' INIT_OPT = 'phi'
init_background = 0 init_background = 1e-05
init_noiselvl = 1e-05 init_noiselvl = 0
iseed = 42 iseed = 42
/ /
&TIME_INTEGRATION_PAR &TIME_INTEGRATION_PAR
......
wk/benchmark_and_scan_scripts/old/lin_KBM.m
+ 133
− 159
View file @ 79443f28
%% QUICK RUN SCRIPT %% QUICK RUN SCRIPT
% This script create a directory in /results and run a simulation directly % This script creates a directory in /results and runs a simulation directly
% from matlab framework. It is meant to run only small problems in linear % from the Matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy" % for benchmarking and debugging purposes since it makes Matlab "busy".
%
% SIMID = 'test_circular_geom'; % Name of the simulation %% Set up the paths for the necessary Matlab modules
% SIMID = 'linear_CBC'; % Name of the simulation gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters
SIMID = 'lin_KBM'; % Name of the simulation SIMID = 'lin_KBM'; % Name of the simulation
RUN = 1 ; % To run or just to load RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/'; EXECNAME = 'gyacomo23_sp'; % single precision
EXECNAME = 'helaz3'; % EXECNAME = 'gyacomo23_dp'; % double precision
% EXECNAME = 'helaz3_dbg';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set up physical parameters
%% Set Up parameters CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss NU = 0.005; % Collision frequency
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TAU = 1.0; % e/i temperature ratio
%% PHYSICAL PARAMETERS K_Ne = 3; % ele Density '''
NU = 0.0; % Collision frequency K_Te = 4.5; % ele Temperature '''
TAU = 1.0; % e/i temperature ratio K_Ni = 3; % ion Density gradient drive
K_Ne = 3; % ele Density ''' K_Ti = 8; % ion Temperature '''
K_Te = 4.5; % ele Temperature '''
K_Ni = 3; % ion Density gradient drive
K_Ti = 8; % ion Temperature '''
SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) 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) NA = 2; % number of kinetic species
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.03; % electron plasma beta BETA = 0.03; % electron plasma beta
%% GRID PARAMETERS %% Set up grid parameters
P = 4; P = 4;
J = P/2; J = P/2;%P/2;
PMAXE = P; % Hermite basis size of electrons PMAX = P; % Hermite basis size
JMAXE = J; % Laguerre " JMAX = J; % Laguerre basis size
PMAXI = P; % " ions NX = 9; % real space x-gridpoints
JMAXI = J; % " NY = 2; % real space y-gridpoints
NX = 9; % real space x-gridpoints LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
NY = 2; % '' y-gridpoints LY = 2*pi/0.25; % Size of the squared frequency domain in y direction
LX = 2*pi/0.1; % Size of the squared frequency domain NZ = 16; % number of perpendicular planes (parallel grid)
LY = 2*pi/0.25; % Size of the squared frequency domain SG = 0; % Staggered z grids option
NZ = 16; % number of perpendicular planes (parallel grid) NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY %% GEOMETRY
% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps % GEOMETRY= 's-alpha';
GEOMETRY= 's-alpha'; GEOMETRY= 'miller';
% GEOMETRY= 'circular'; EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear SHEAR = 0.8; % magnetic shear
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) KAPPA = 1.0; % elongation
EPS = 0.18; % inverse aspect ratio DELTA = 0.0; % triangularity
%% TIME PARMETERS ZETA = 0.0; % squareness
TMAX = 15.1; % Maximal time unit PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
DT = 5e-3; % Time step SHIFT_Y = 0.0; % Shift in the periodic BC in z
SPS0D = 1; % Sampling per time unit for 2D arrays NPOL = 1; % Number of poloidal turns
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 10; % Sampling per time unit for 2D arrays %% TIME PARAMETERS
SPS5D = 1/5; % Sampling per time unit for 5D arrays TMAX = 15; % Maximal time unit
SPSCP = 0; % Sampling per time unit for checkpoints DT = 5e-3; % Time step
JOB2LOAD= -1; 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 %% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) GKCO = 0; % Gyrokinetic operator
CO = 'DG'; ABCO = 1; % INTERSPECIES collisions
GKCO = 0; % gyrokinetic operator INIT_ZF = 0; % Initialize zero-field quantities
ABCO = 1; % INTERSPECIES collisions HRCY_CLOS = 'truncation'; % Closure model for higher order moments
INIT_ZF = 0; ZF_AMP = 0.0; DMAX = -1;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) NMAX = 0;
KERN = 0; % Kernel model (0 : GK) KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
%% OUTPUTS %% OUTPUTS
W_DOUBLE = 1; W_DOUBLE = 1; % Output flag for double moments
W_GAMMA = 1; W_HF = 1; W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
W_PHI = 1; W_NA00 = 1; W_HF = 1; % Output flag for high-frequency potential energy
W_DENS = 1; W_TEMP = 1; W_PHI = 1; % Output flag for potential
W_NAPJ = 1; W_SAPJ = 0; W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
W_DENS = 1; % Output flag for total density
W_TEMP = 1; % Output flag for temperature
W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused %% UNUSED PARAMETERS
HD_CO = 0.0; % Hyper diffusivity cutoff ratio % These parameters are usually not to play with in linear runs
MU = 0.0; % Hyperdiffusivity coefficient MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % Hyperdiffusivity coefficient in x direction
MU_X = MU; % MU_Y = MU; % Hyperdiffusivity coefficient in y direction
MU_Y = MU; % N_HD = 4; % Degree of spatial-hyperdiffusivity
N_HD = 4; MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
MU_Z = 10; % HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
MU_P = 0.0; % MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
MU_J = 0.0; % MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
LAMBDAD = 0.0; LAMBDAD = 0.0; % Lambda Debye
NOISE0 = 0.0e-5; % Init noise amplitude NOISE0 = 0.0e-5; % Initial noise amplitude
BCKGD0 = 1.0; % Init background BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator
%%------------------------------------------------------------------------- %%-------------------------------------------------------------------------
%% RUN %% RUN
setup setup
% system(['rm fort*.90']); % system(['rm fort*.90']);
% Run linear simulation % Run linear simulation
if RUN if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) % RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end end
%% Load results %% Analysis
%% % load
filename = [SIMID,'/',PARAMS,'/']; filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
LOCALDIR = [HELAZDIR,'results/',filename,'/']; LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
% Load outputs from jobnummin up to jobnummax % Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00; J0 = 0; J1 = 0;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
%% Short analysis
if 0 if 0
%% linear growth rate (adapted for 2D zpinch and fluxtube) %% Plot heat flux evolution
trange = [0.5 1]*data.Ts3D(end); figure
nplots = 3; semilogy(data.Ts0D,data.HFLUX_X);
lg = compute_fluxtube_growth_rate(data,trange,nplots); xlabel('$tc_s/R$'); ylabel('$Q_x$');
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
end
if 1
%% Ballooning plot
options.time_2_plot = [120];
options.kymodes = [0.25];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
end end
if 1 % Activate or not
if 0 %% plot mode evolution and growth rates
%% Hermite-Laguerre spectrum % Load phi
% options.TIME = 'avg'; [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
options.P2J = 1; options.NORMALIZED = 0;
options.ST = 1; options.TIME = data.Ts3D;
options.PLOT_TYPE = 'space-time'; % Time window to measure the growth of kx/ky modes
% options.PLOT_TYPE = 'Tavg-1D'; options.KX_TW = [0.2 1]*data.Ts3D(end);
% options.PLOT_TYPE = 'Tavg-2D'; options.KY_TW = [0.2 1]*data.Ts3D(end);
options.NORMALIZED = 0; options.NMA = 1; % Set NMA option to 1
options.JOBNUM = 0; options.NMODES = 999; % Set how much modes we study
options.TIME = [0 50]; options.iz = 'avg'; % Compressing z
options.specie = 'i'; options.ik = 1; %
options.compz = 'avg'; options.fftz.flag = 0; % Set fftz.flag option to 0
fig = show_moments_spectrum(data,options); fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
% fig = show_napjz(data,options);
save_figure(data,fig)
end end
if 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
options.kzkx = 0;
options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = [0:1000];
options.NMA = 1;
options.NMODES = 1;
options.iz = 'avg';
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
if 0
%% RH TEST
ikx = 2; t0 = 0; t1 = data.Ts3D(end);
[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3));
figure
plot(data.Ts3D(it0:it1), plt(data.PHI));
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
end
wk/benchmark_and_scan_scripts/plot_sparc_edge_profiles.m 0 → 100644
+ 30
− 0
View file @ 79443f28
d_ = load("rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt");
rho_ne = d_(:,1); ne = d_(:,2);
d_ = load("rodriguez_fernandez_2020_fig3_sparc_T_prof.txt");
rho_Te = d_(:,1); Te = d_(:,2);
figure;
plot(rho_ne,ne,'.k','DisplayName','$n_e$'); hold on;
plot(rho_Te,Te,'or','DisplayName','$T_e$');
%%
[~,i1] = min(abs(rho_ne-0.95));
[~,i2] = min(abs(rho_ne-0.98));
[~,in975] = min(abs(rho_ne-0.975));
fit_n = polyfit(rho_ne(i1:i2),ne(i1:i2),1);
rho_ = 0.925:0.025:1.0;
plot(rho_,fit_n(1).*rho_+fit_n(2),'-k','DisplayName',['$-\partial_x n/n_{975}=$',num2str(-fit_n(1)/ne(in975))]);
%%
[~,i1] = min(abs(rho_Te-0.95));
[~,i2] = min(abs(rho_Te-0.98));
[~,iT975] = min(abs(rho_Te-0.975));
fit_T = polyfit(rho_Te(i1:i2),Te(i1:i2),1);
rho_ = 0.925:0.025:1.0;
plot(rho_,fit_T(1).*rho_+fit_T(2),'-r','DisplayName',['$-\partial_x T/T_{975}=$',num2str(-fit_T(1)/Te(iT975))]);
xline(rho_Te(iT975));
legend('show')
\ No newline at end of file
wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_T_prof.txt 0 → 100644
+ 68
− 0
View file @ 79443f28
0.7368446165421814, 6.290726909030852
0.7448753005083408, 6.144500926758546
0.7468572791837185, 6.364353746490245
0.75290231414362, 6.034904846671925
0.758917986456479, 5.998495164339065
0.768949000752418, 5.888972490870049
0.7809840157090162, 5.779523224018643
0.7870033583527556, 5.7064836395000995
0.7950046796718726, 5.853296874713255
0.795030371988035, 5.5968875594134815
0.8030353636380321, 5.707070892440953
0.8110880695894738, 5.3410654970545615
0.8151070819034336, 5.231322603732725
0.8250977225596889, 5.524728854306215
0.8271237452056306, 5.305022847809731
0.8351580995026704, 5.122166963351745
0.837140078178048, 5.3420197830834475
0.8411774421464093, 5.049127378833202
0.8532197977647686, 4.866418307610431
0.8552164577636674, 4.939751518599397
0.8592317997467473, 4.866638527463252
0.8632398010680664, 4.866785340698463
0.8692554733809255, 4.830375658365604
0.8913251729643428, 4.574773815859498
0.9073681892422603, 4.465471362243306
0.9113761905635795, 4.465618175478522
0.9153731808922576, 4.575654695270778
0.915387862215779, 4.429135086528049
0.9294342184947975, 4.246499421922888
0.9294342184947975, 4.246499421922888
0.9354425501458956, 4.283349543961386
0.945477234772715, 4.137196968306689
0.9494925767557947, 4.064083977170544
0.953533611055037, 3.7345616707346174
0.957508579398433, 4.0643776036409704
0.957508579398433, 4.0643776036409704
0.957552623368997, 3.6248187774127842
0.9595309317134941, 3.881301499330167
0.9595676350222972, 3.5150024774733453
0.9596043383311006, 3.1487034556165234
0.9596116789928613, 3.075443651245159
0.9615642950211962, 3.5883356884623154
0.961652382962324, 2.7092180360059466
0.9635866473362575, 3.405259584151512
0.9636784056082659, 2.4895120295094593
0.9656163403130793, 3.1489236754693444
0.9656493732910023, 2.8192545557982065
0.9656970875924467, 2.34306582738434
0.9677194399075078, 2.1599897230735365
0.969686737259364, 2.526362151547964
0.9697528032152098, 1.8670239122056849
0.97165036428034, 2.9293644822080687
0.9737094199042047, 2.379989356040447
0.9737975078453325, 1.5008717035840782
0.9757354425501461, 2.1602833495439633
0.9758125194986329, 1.3910554036446392
0.977757794865207, 1.9772072452331564
0.9777724761887283, 1.8306876364904312
0.9818024994953296, 1.6110550366115497
0.9818355324732526, 1.2813859169404118
0.9818648951202953, 0.9883466994549543
0.9838578847883137, 1.098309812629605
0.9879062597493165, 0.6955277018223178
0.9879062597493165, 0.6955277018223178
0.9898882384246943, 0.9153805215540203
0.989895579086455, 0.8421207171826559
0.9939182617312952, 0.6957479216751388
0.9999449450367952, 0.549448532785231
wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt 0 → 100644
+ 43
− 0
View file @ 79443f28
0.8907345182249297, 2.9380728997180983
0.9087345840386986, 2.9383361547939852
0.9147309496561232, 2.9237986244922283
0.9227236662656902, 2.894665062760742
0.9307236955162542, 2.8947820650166918
0.9367200611336789, 2.8802445347149357
0.9407164194384622, 2.8656777538491927
0.9447091214227475, 2.83648569098973
0.9486945107660357, 2.778043064142831
0.9526616185068318, 2.64647402732734
0.9526799001093241, 2.7196004372959317
0.9526799001093241, 2.7196004372959317
0.95664700785012, 2.588031400480441
0.9585994829962815, 2.397931985126089
0.9586323898807673, 2.5295595230695547
0.9605885213474273, 2.3540853897089207
0.9625592780960809, 2.2371123843231615
0.962570247057576, 2.280988230304315
0.9625848723395699, 2.33948935827919
0.9645263785242358, 2.1055140969436827
0.9645300348447342, 2.1201393789374023
0.9645519727677249, 2.207891070899713
0.9645519727677249, 2.207891070899713
0.9665117605548833, 2.047042219532795
0.9684825173035368, 1.9300692141470357
0.9685007989060289, 2.0031956241156283
0.9704569303726888, 1.8277214907549926
0.9723509043908755, 1.4036175635011467
0.9723655296728692, 1.4621186914760216
0.9723801549548629, 1.5206198194508946
0.9723911239163582, 1.5644956654320499
0.9724057491983519, 1.6229967934069247
0.9724203744803456, 1.6814979213817978
0.9724386560828379, 1.7546243313503904
0.9763106994906746, 1.2427979626982193
0.9763289810931667, 1.3159243726668102
0.9802741509109724, 1.0966036438890097
0.9802851198724677, 1.140479489870165
0.9822558766211212, 1.023506484484404
0.9842485712927652, 0.9942851710609553
0.9842485712927652, 0.9942851710609553
0.9882266479950569, 0.9065919802266187
0.9882303043155553, 0.9212172622203383
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