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

matlab scripts update

parent dc225672
No related branches found
No related tags found
No related merge requests found
Showing
with 249 additions and 249 deletions
% HermitePoly.m by David Terr, Raytheon, 5-10-04
% Given nonnegative integer n, compute the
% Hermite polynomial H_n. Return the result as a vector whose mth
% element is the coefficient of x^(n+1-m).
% polyval(HermitePoly(n),x) evaluates H_n(x).
function hk = HermitePoly(n)
% Evaluate the normalized Hermite polynomial.
if n==0
hk = 1;
elseif n==1
hk = [2 0];
else
hkm2 = zeros(1,n+1);
hkm2(n+1) = 1;
hkm1 = zeros(1,n+1);
hkm1(n) = 2;
for k=2:n
hk = zeros(1,n+1);
for e=n-k+1:2:n
hk(e) = 2*(hkm1(e+1) - (k-1)*hkm2(e));
end
hk(n+1) = -2*(k-1)*hkm2(n+1);
if k<n
hkm2 = hkm1;
hkm1 = hk;
end
end
end
\ No newline at end of file
function [SS,XX,FF] = compute_fa_2D(data, species, s, x, T)
function [SS,XX,FF,FAM] = compute_fa_2D(data, species, s, x, T)
%% Compute the dispersion function from the moment hierarchi decomp.
% Normalized Hermite
Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p));
% Hp = @(p,s) hermiteH(p,s)./sqrt(2.^p.*factorial(p));
% Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p));
Hp = @(p,s) polyval(HermitePhys_norm(p),s);
% Laguerre
Lj = @(j,x) polyval(LaguerrePoly(j),x);
% Maxwellian factor
......@@ -18,30 +18,22 @@ switch species
case 'i'
Napj_ = data.Napjz(1,:,:,:,:);
end
parray = double(data.grids.Parray);
jarray = double(data.grids.Jarray);
% switch options.iz
% case 'avg'
options.SHOW_FLUXSURF = 0;
options.SHOW_METRICS = 0;
options.SHOW_CURVOP = 0;
[~, geo_arrays] = plot_metric(data,options);
J = geo_arrays.Jacobian;
Nz = data.grids.Nz;
tmp_ = 0;
for iz = 1:Nz
tmp_ = tmp_ + J(iz)*Napj_(:,:,:,iz,:);
end
Napj_ = tmp_/sum(J(iz));
% Napj_ = mean(Napj_,4);
% Napj_ = Napj_(:,:,:,Nz/2+1,:);
% phi_ = mean(data.PHI,3);
% otherwise
% iz = options.iz;
% Napj_ = Napj_(:,:,:,:,iz,:);
% phi_ = data.PHI(:,:,iz);
% end
% Napj_ = squeeze(Napj_);
options.SHOW_FLUXSURF = 0;
options.SHOW_METRICS = 0;
options.SHOW_CURVOP = 0;
[~, geo_arrays] = plot_metric(data,options);
J = geo_arrays.Jacobian;
Nz = data.grids.Nz;
tmp_ = 0;
for iz = 1:Nz
tmp_ = tmp_ + J(iz)*Napj_(:,:,:,iz,:);
end
Napj_ = tmp_/sum(J(iz));
frames = T;
for it = 1:numel(T)
......@@ -56,43 +48,19 @@ Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
% if options.RMS
FF = zeros(numel(x),numel(s));
FAM = FaM(SS,XX);
for ip_ = 1:Np
p_ = parray(ip_);
HH = Hp(p_,SS);
for ij_ = 1:Nj
j_ = jarray(ij_);
LL = Lj(j_,XX);
HLF = HH.*LL.*FAM;
FF = FF + Napj_(ip_,ij_)*HLF;
end
end
% else
% FF = zeros(numel(options.XPERP),numel(options.SPAR));
% FAM = FaM(SS,XX);
% for ip_ = 1:Np
% p_ = parray(ip_);
% HH = Hp(p_,SS);
% for ij_ = 1:Nj
% j_ = jarray(ij_);
% LL = Lj(j_,XX);
% FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM;
% end
% end
% end
FF = zeros(numel(x),numel(s));
FAM = FaM(SS,XX);
for ip_ = 1:Np
p_ = parray(ip_);
HH = Hp(p_,SS);
for ij_ = 1:Nj
j_ = jarray(ij_);
LL = Lj(j_,XX);
HLF = HH.*LL.*FAM;
FF = FF + Napj_(ip_,ij_)*HLF;
end
end
FF = (FF.*conj(FF)); %|f_a|^2
% FF = abs(FF); %|f_a|
% if options.RMS
% FF = squeeze(mean(mean(sqrt(FF),1),2)); %sqrt(<|f_a|^2>kx,ky)
FF = sqrt(FF); %<|f_a|>kx,ky
% else
% FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y)
% end
% FF = FF./max(max(FF));
% FF = FF';
% FF = sqrt(FF);
% FF = FF';
FF = sqrt(FF); %<|f_a|>kx,ky
FF = FF./max(FF(:));
end
\ No newline at end of file
......@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
tw = [00 1000];
tw = [100 1000];
fig = gcf;
axObjs = fig.Children;
......@@ -29,8 +29,11 @@ figure;
plot(mvm(X_(n0:skip:n1)-shift),mvm(Y_(n0:skip:n1))); hold on;
% t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
avg= mean(Y_(n0:n1)); dev = std(Y_(n0:n1));
[avg, sliceav, sliceerr] = sliceAverage(mvm(Y_(n0:skip:n1)),10);
dev = std(sliceav);
% avg= mean(Y_(n0:n1)); dev = std(Y_(n0:n1));
disp(['AVG =',sprintf('%4.4f',avg),'+-',sprintf('%4.4f',dev)]);
disp([sprintf('%4.4f',avg),',',sprintf('%4.4f',dev)]);
%
% n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
% avg_ = mean(Y_(n1:n2));
......
......@@ -104,11 +104,5 @@ for i_ = 1:numel(names)
namae = names{i_};
arrays.(namae) = geo_arrays(:,i_);
end
try
if ~options.SHOWFIG
close
end
catch
end
end
function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
function [ FIGURE, XX, YY, ZZ ] = show_moments_spectrum( DATA, OPTIONS )
species_name = {'i','e'}; % hard coded because this list is not output yet
XX = 0; YY = 0; ZZ = 0;
Pa = DATA.grids.Parray;
Ja = DATA.grids.Jarray;
Time_ = DATA.Ts3D;
......@@ -67,6 +67,9 @@ for ia = 1:DATA.inputs.Na
[~,it0] = min(abs(DATA.Ts3D-t0));
[~,it1] = min(abs(DATA.Ts3D-t1));
Napjz_tavg = mean(Napjz(:,:,it0:it1),3);
if OPTIONS.NORMALIZED
Napjz_tavg = Napjz_tavg./max(Napjz_tavg(:));
end
x_ = DATA.grids.Parray;
y_ = DATA.grids.Jarray;
if OPTIONS.TAVG_2D_CTR
......@@ -77,7 +80,7 @@ for ia = 1:DATA.inputs.Na
set(gca,'YDir','normal')
xlabel('$p$'); ylabel('$j$');
end
clb = colorbar; colormap(jet);
clb = colorbar; colormap(parula);
clb.Label.String = '$\langle E(p,j) \rangle_t$';
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
......
......@@ -189,6 +189,7 @@ if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false
if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end;
if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end;
if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end;
if W_FVEL; OUTPUTS.write_fvel = '.true.'; else; OUTPUTS.write_fvel = '.false.';end;
if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.false.';end;
OUTPUTS.job2load = JOB2LOAD;
%% Create directories
......
......@@ -54,6 +54,7 @@ fprintf(fid,[' write_Na00 = ', OUTPUTS.write_Na00,'\n']);
fprintf(fid,[' write_Napj = ', OUTPUTS.write_Napj,'\n']);
fprintf(fid,[' write_dens = ', OUTPUTS.write_dens,'\n']);
fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']);
fprintf(fid,[' write_fvel = ', OUTPUTS.write_fvel,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&MODEL\n');
......
PARTITION = '/misc/gyacomo23_outputs/paper_3/';
switch 3
PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
switch 4
case 1
SIM_SET_NAME = 'Multi-scale';
E_FLUX = 1;
FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_3x2x768x192x24';
% FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_5x2x768x192x24';
case 2
SIM_SET_NAME = 'Ion-scale';
SIM_SET_NAME = 'KEM';
E_FLUX = 1;
FOLDER = 'DIIID_fullphys_rho95_geom_scan/ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0';
FOLDER = 'ion_scale/5x2x256x64x32';
case 3
SIM_SET_NAME = 'Adiab. e.';
SIM_SET_NAME = 'AEM';
E_FLUX = 0;
FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/';
FOLDER = 'adiabatic_electrons/5x2x256x64x32';
case 4
SIM_SET_NAME = 'Adiab. e.';
SIM_SET_NAME = 'RFM';
E_FLUX = 0;
FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/';
FOLDER = 'hot_electrons/256x64x32';
end
GEOM = {'NT','0T','PT'};
......
......@@ -4,72 +4,9 @@ addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/rm_corrections_HF/';
% folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2a_gyroLES/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_5e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_gyroLES/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_1e-2_muv_1e-1/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
% folder = '/misc/gene_results/Z-pinch/kN_2.0_HD_transport_spectrum_01/';
% folder = '/misc/gene_results/Z-pinch/kN_2.5_HD_transport_spectrum_01/';
folder = '/misc/gyacomo23_outputs/triangularity_paper/GENE_baseline/output/';
% folder = '/misc/gyacomo23_outputs/triangularity_paper/GENE_output/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
% folder = '/misc/gene_results/CBC/256x96x24x32x12/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
% folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
% folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
% folder = '/misc/gene_results/miller/';
% folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/KT_5.3_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
% folder = '/misc/gene_results/linear_miller_CBC/hdf5_miller_s0_adiabe/';
% folder = '/misc/gene_results/linear_miller_CBC/hdf5_salpha_s0_adiabe/';
%Paper 2
% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
% folder = '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_4.0/';
% folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
% Miller CBC
% folder = '/misc/gene_results/Miller_AUG/Npol_1/output/';
% folder = '/misc/gene_results/Miller_AUG/Npol_5/output/';
% folder = '/misc/gene_results/Miller_AUG/circ/';
folder = '/misc/gene_results/Miller_AUG/linear/Npol_1/';
% debug ? shearless
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_00/';
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_01/';
if 1
%% FULL DATA LOAD (LONG)
gene_data = load_gene_data(folder);
......@@ -163,14 +100,14 @@ options.NAME = '\phi';
options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.COMP = 17;
options.TIME = gene_data.Ts3D;
options.RESOLUTION= 256;
gene_data.a = data.EPS * 2000;
% gene_data.a = data.EPS * 2000;
create_film(gene_data,options,'.gif')
end
if 1
if 0
%% Geometry
names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
'$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
......
......@@ -5,28 +5,48 @@ addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add
default_plots_options
% Partition of the computer where the data have to be searched
% PARTITION='/Users/ahoffmann/gyacomo/results/paper_3/';
% PARTITION='/Users/ahoffmann/gyacomo/results/';
PARTITION='/home/ahoffman/gyacomo/results/';
% PARTITION='/misc/gyacomo23_outputs/paper_3/';
% resdir = 'DIIID_rho_95_Tstudy/ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/0T/';
% resdir = 'DIIID_rho_95_Tstudy/adiab_e/5x2x256x64x32_tau_1_RN_0/0T/';
%%
% resdir = 'AE_3x2x128x32x24/PT';
% resdir = 'AE_3x2x128x32x24/PT';
% resdir = 'AE_5x3x128x32x24/NT';
resdir = 'IS_5x3x128x32x24/PT';
% Triangularity paper
PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
% resdir = 'ion_scale/3x2x256x64x32/PT';
% resdir = 'ion_scale/5x2x256x64x32/0T';
% Nominal parameters
% resdir = 'ion_scale/3x2x256x64x32/0T';
% resdir = 'ion_scale/5x3x256x64x32/0T';
% resdir = 'ion_scale/5x3x192x48x24/0T';
resdir = 'ion_scale/9x5x256x64x32/0T';
% resdir = 'ion_scale/restart/5x3x256x64x32/0T';
% resdir = 'ion_scale/restart/9x5x192x48x24/0T';
% resdir = 'adiabatic_electrons/5x2x256x64x32/0T';
% resdir = 'adiabatic_electrons/5x2x192x48x24/NT';
% resdir = 'adiabatic_electrons/5x2x192x48x24/0T';
% resdir = 'hot_electrons/256x64x32/0T';
% resdir = 'hot_electrons/L_300/256x64x32/0T';
resdir = 'hot_electrons/L_300_gradN_scaled/256x64x32/0T';
% resdir = 'hot_electrons/256x64x32/0T';
% resdir = 'hot_electrons/512x64x32/0T';
% PARTITION = '/misc/gyacomo23_outputs/paper_3/';
% resdir = 'DIIID_HEL_rho95/PT';
% No grad N
% PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/no_gradN';
% resdir = '/ion_scale/3x2x256x64x24/0T';
% resdir = '/ion_scale/5x2x256x64x24/0T';
% resdir = '/ion_scale/9x5x128x32x24/0T';
% resdir = '/adiabatic_electrons/3x2x256x64x24/0T';
% resdir = '/adiabatic_electrons/5x2x256x64x24/0T';
% resdir = '/hot_electrons/L_300/256x64x32/0T';
% PARTITION = '/home/ahoffman/gyacomo/results/lin_DIIID_LM_rho95_scan/';
% resdir = '6x2x32_5x3_Lx_120_Ly_37.6991_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
% resdir = '6x2x32_17x9_Lx_120_Ly_37.6991_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
% resdir = '6x2x32_5x3_Lx_120_Ly_12.5664_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
% resdir = '6x2x32_17x9_Lx_120_Ly_12.5664_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
% resdir = '6x2x32_5x3_Lx_120_Ly_8.1955_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
% resdir = '6x2x32_17x9_Lx_120_Ly_8.1955_q0_4.8_e_0.3_s_2.5_mil__kN_1.7_kT_5.2_nu_2.0e-02_DGGK/';
DATADIR = [PARTITION,resdir,'/'];
% read_flux_out_XX(DATADIR,1,5)
read_flux_out_XX(DATADIR,1,1);
%%
J0 = 00; J1 = 10;
......@@ -55,7 +75,7 @@ if data.inputs.Na > 1
data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
end
end
if 1
if 0
%% Plot transport and phi radial profile
% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
......@@ -79,7 +99,7 @@ options.RESOLUTION = 256;
plot_radial_transport_and_spacetime(data,options);
end
if 1
if 0
%% 2D field snapshots
% Options
options.INTERP = 0;
......@@ -89,13 +109,13 @@ options.NORMALIZE = 0;
options.LOGSCALE = 0;
options.CLIMAUTO = 1;
options.TAVG = 1;
% options.NAME = ['N_i^{00}'];
options.NAME = ['N_i^{00}'];
% options.NAME = 'n_e';
% options.NAME = 'u_i';
% options.NAME = 'n_i';
% options.NAME = 'Q_{xi}';
% options.NAME = 'v_{Ey}';
options.NAME = 'w_{Ez}';
% options.NAME = 'w_{Ez}';
% options.NAME = '\omega_z';
% options.NAME = '\phi';
% options.NAME = 'n_i-n_e';
......@@ -103,12 +123,12 @@ loc =11;
[~,i_] = min(abs(loc - data.grids.y));
options.COMP =i_;
% options.PLAN = '3D';
options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
% options.PLAN = 'xz'; options.COMP ='avg';
% options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
options.PLAN = 'xz'; options.COMP ='avg';
% options.COMP ='avg';
options.XYZ =[-11 20 0];
options.TIME = [100 250 350 500]; options.TAVG = 0;
% options.TIME = [100:250]; options.TAVG = 1;
% options.TIME = [100 250 350 500]; options.TAVG = 0;
options.TIME = [100:500]; options.TAVG = 1;
options.RESOLUTION = 256;
fig = photomaton(data,options);
% colormap(gray)
......@@ -120,7 +140,7 @@ if 0
profiler(data)
end
if 1
if 0
%% Mode evolution
% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
......@@ -129,8 +149,8 @@ if 1
options.NORMALIZED = 0;
options.TIME = data.Ts3D;
options.KX_TW = [0.1 2.5]; %kx Growth rate time window
options.KY_TW = [0.1 2.5]; %ky Growth rate time window
options.KX_TW = [0.5 1]*data.Ts3D(end); %kx Growth rate time window
options.KY_TW = [0.5 1]*data.Ts3D(end); %ky Growth rate time window
options.NMA = 1;
options.NMODES = 64;
options.iz = 'avg'; % avg or index
......@@ -173,19 +193,19 @@ options.SHOWFIG = 1;
% % save_figure(data,fig,'.png')
end
if 1
if 0
%% Hermite-Laguerre spectrum
[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
data.Napjz(1,3,1,:,:) = data.Napjz(1,3,1,:,:)*data.inputs.tau;
data.Napjz(1,1,2,:,:) = data.Napjz(1,1,2,:,:)*data.inputs.tau;
data.Napjz(1,3,1,:,:) = data.Napjz(1,3,1,:,:)*data.inputs.tau(1);
data.Napjz(1,1,2,:,:) = data.Napjz(1,1,2,:,:)*data.inputs.tau(1);
% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 0;
options.NORMALIZED = 1;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
options.TAVG_2D_CTR= 0; %make it contour plot
options.TWINDOW = [6 20];
options.TWINDOW = [20 20];
fig = show_moments_spectrum(data,options);
end
......@@ -204,15 +224,15 @@ options.SHOWFIG = 1;
[fig, chi, phib, psib, ~] = plot_ballooning(data,options);
end
if 1
if 0
%% 1D spectral plot
options.TIME = [100 300]; % averaging time window
% options.NAME = ['N_i^{00}'];
options.NAME = ['N_i^{00}'];
% options.NAME = 'n_i';
% options.NAME = 'T_i';
% options.NAME = 'Q_{xi}';
% options.NAME = 's_{Ey}';
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = '\psi';
options.NORMALIZE = 0;
[fig] = plot_spectrum(data,options);
......@@ -231,7 +251,7 @@ if 0
% data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
% [data.TEMP, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'temp');
% data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 0;
options.INTERP = 1;
options.POLARPLOT = 0;
options.BWR = 1; % bluewhitered plot or gray
options.CLIMAUTO = 0; % adjust the colormap auto
......@@ -240,8 +260,8 @@ options.CLIMAUTO = 0; % adjust the colormap auto
% options.NAME = '\psi';
% options.NAME = 'n_i';
% options.NAME = '\phi^{NZ}';
% options.NAME = ['N_e^{00}'];
options.NAME = ['N_i^{00}'];
% options.NAME = ['N_i^{00}'];
options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
% options.PLAN = 'xz'; options.COMP ='avg';
% options.PLAN = '3D';
......@@ -275,4 +295,10 @@ for i = 1:nSV
'color',colors_(i,:),'DisplayName',['SV ',num2str(i)]);hold on
end
legend('show');
end
if 0
%% Pseudo fluid analysis
Time_window = [100 200];
pseudo_fluid_analysis
end
\ No newline at end of file
......@@ -48,20 +48,26 @@ end
% SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)];
%% Change parameters
% GEOMETRY = 's-alpha';
PMAX = 4; JMAX = 2;
DELTA = 0.2;
% K_tB = 0; K_mB = 0; K_ldB = 0;
% K_Ni = 0; K_Ne = 0;
% K_Ti = 0; TAU = 1/3;
PMAX = 8; JMAX = 4;
% DELTA = 0.2;
% K_tB = 0; K_mB = 0; K_ldB = 1;
% K_Ni = 0;
% K_Ne = 0;
% K_Ti = 0;
% K_Te = 0;
DELTA = 0.0;
% DELTA = 0.2;
S_DELTA = DELTA/2;
LY = 2*pi/0.75;
LY = 2*pi/0.1;
TMAX = 40;
NY = 2;
DT = 0.0025;
NY = 20;
DT = 0.001;
CO = 'DG';
% TAU = 1; NU = 0.05;
% TAU = 1e-3; K_Ti = K_Ti/2/TAU; NU = 3*NU/8/TAU; ADIAB_E = 1; NA = 1;
NA = 1; ADIAB_E = 1; DT = 1e-2; DTSAVE3D = 1e-2; TMAX = 60;
TAU = 1e-3; K_Ti = K_Ti/2/TAU; K_Ni = 0;
NU = 3*NU/8/TAU; PMAX = 2; JMAX = 1;
% K_Ti = K_Ti/4;
% MU_X = 1; MU_Y = 1;
%% RUN
setup
......@@ -201,4 +207,19 @@ if 0
subplot(1,2,2); plot(R,Z,'-k'); xlabel('$R$'); ylabel('$Z$');
axis equal
title(data.paramshort);
end
\ No newline at end of file
end
if 0
%% Hermite-Laguerre spectrum
[data.Napjz, data.Ts3D] = compile_results_3Da(LOCALDIR,J0,J1,'Napjz');
% data.Napjz(1,3,1,:,:) = data.Napjz(1,3,1,:,:)*data.inputs.tau;
% data.Napjz(1,1,2,:,:) = data.Napjz(1,1,2,:,:)*data.inputs.tau;
% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
options.ST = 1;
options.NORMALIZED = 1;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
options.TAVG_2D_CTR= 0; %make it contour plot
options.TWINDOW = [20 20];
fig = show_moments_spectrum(data,options);
end
......@@ -39,7 +39,7 @@ run lin_DIIID_LM_rho95
% NU = 1;
% TAU = 1;
NY = 2;
DELTA =-0.2; TRIANG = 'NT';
DELTA =0.0; TRIANG = '';
S_DELTA = DELTA/2;
% EXBRATE = 0;
% S_DELTA = min(2.0,S_DELTA);
......@@ -48,15 +48,23 @@ S_DELTA = DELTA/2;
LX = 120;
%% Scan parameters
SIMID = [SIMID,TRIANG,'_scan'];
P_a = [2 4 6 8 16]; J_a = [1 2 3 4 8];
P_a = [2 4 8 16]; J_a = [1 2 4 8];
% P_a = [2 4]; J_a = [1 1];
% P_a = 2;
% ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0];
ky_a = [0.05 linspace(0.1,1.1,16)]; ky_a = ky_a(1:end-2);
% ky_a = 4.0;
% dt_a = logspace(-2,-3,numel(ky_a));
DT = 0.0005;
CO = 'DG';
DTSAVE3D = 0.002; TMAX = 40;
% KEM
NA = 2; ADIAB_E = 0; DT = 5e-4; DTSAVE3D = 5e-3; TMAX = 60;
% AEM
% NA = 1; ADIAB_E = 1; DT = 1e-3; DTSAVE3D = 5e-2; TMAX = 60;
%RFM
% NA = 1; ADIAB_E = 1; DT = 5e-3; DTSAVE3D = 1e-2; TMAX = 60;
% TAU = 1e-3; K_Ti = K_Ti/2/TAU; K_Ni = 0;
% NU = 3*NU/8/TAU; P_a = 2; J_a = 1; ky_a = 2*ky_a;
K_Ni = 0;
%% Scan loop
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a));
......
......@@ -10,13 +10,27 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.32141_LDGK_0.0038027_be_0.0060039.mat';
% datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.0883_LDGK_0.0080915_be_0.0015991.mat';
% rho = 0.95
% datagname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.62888_LDGK_0.0046858_be_0.0048708.mat';
% datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_4_kN_0.62888_DGGK_0.0046858_be_0.0048708.mat';
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.025_1.1_P_2_4_kN_1.7_DGGK_0.02_be_0.000759_d_0.0.mat';
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.025_1.1_P_2_6_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.1_1.1_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
% KEM
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
datafname = 'lin_DIIID_LM_rho95PT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.2.mat';
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.05_0.96667_P_2_4_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
% datafname = 'lin_DIIID_LM_rho95PT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.2.mat';
% datafname = 'lin_DIIID_LM_rho95NT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_-0.2.mat';
% AEM
% datafname = 'lin_DIIID_LM_rho95AE_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
% datafname = 'lin_DIIID_LM_rho95PTAE_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.2.mat';
% datafname = 'lin_DIIID_LM_rho95NTAE_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_-0.2.mat';
% RFM
% datafname = 'lin_DIIID_LM_rho95RFM_scan/6x32_ky_0.1_1.9333_P_2_2_kN_0_DGGK_7.5_be_0.000759_d_0.mat';
%
% no gradn
% KEM
% datafname = 'lin_DIIID_LM_rho95NT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_0_DGGK_0.02_be_0.000759_d_-0.2.mat';
% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.05_0.96667_P_2_16_kN_0_DGGK_0.02_be_0.000759_d_0.mat';
% datafname = 'lin_DIIID_LM_rho95PT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_0_DGGK_0.02_be_0.000759_d_0.2.mat';
% AEM
% datafname = 'lin_DIIID_LM_rho95AE_scan/6x32_ky_0.05_0.96667_P_2_16_kN_0_DGGK_0.02_be_0.000759_d_0.mat';
% RFM
% datafname = '';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
......@@ -26,7 +40,7 @@ d = load(fname);
gamma = real(d.data); g_err = real(d.err);
omega = imag(d.data); w_err = imag(d.err);
if FILTERGAMMA
gamma = gamma.*(gamma>0.025);
gamma = gamma.*(gamma>0.005);
end
if 0
%% Pcolor of the peak
......@@ -53,8 +67,10 @@ figure
colors_ = jet(numel(d.s2));
subplot(121)
for i = 1:numel(d.s2)
% toplot = ~isnan(gamma(:,i));
toplot = (gamma(:,i)-g_err(:,i)>0);
% plot(d.s1,gamma(:,i),'s-',...
plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
plot(d.s1(toplot),gamma(toplot,i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
......@@ -69,7 +85,9 @@ xlim([d.s1(1) d.s1(end)]);
subplot(122)
for i = 1:numel(d.s2)
plot(d.s1,omega(:,i),'s-',...
toplot = ~isnan(gamma(:,i));
toplot = (gamma(:,i)-g_err(:,i)>0);
plot(d.s1(toplot),omega(toplot,i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
......
% Nominal parameters
PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
% resdir = 'ion_scale/3x2x256x64x32/0T';
% resdir = 'ion_scale/5x2x256x64x32/NT';
% resdir = 'ion_scale/5x3x256x64x32/NT';
% resdir = 'ion_scale/5x3x192x48x24/NT';
resdir = 'ion_scale/9x5x256x64x32/0T';
% resdir = 'ion_scale/restart/5x3x256x64x32/0T';
% resdir = 'ion_scale/restart/9x5x192x48x24';
% resdir = 'adiabatic_electrons/5x2x256x64x32/0T';
% resdir = 'adiabatic_electrons/5x2x192x48x24/0T';
% resdir = 'hot_electrons/256x64x32/0T';
% resdir = 'hot_electrons/256x64x32/0T';
% resdir = 'hot_electrons/512x64x32/0T';
DATADIR = [PARTITION,resdir,'/'];
read_flux_out_XX(DATADIR,1,1);
%%
% No grad N
PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/no_gradN/';
% resdir = '/ion_scale/3x2x256x64x24/0T';
% resdir = '/ion_scale/5x2x256x64x24/PT';
% resdir = '/ion_scale/5x3x256x64x24/PT';
% resdir = '/ion_scale/5x3x192x48x24/NT';
% resdir = '/adiabatic_electrons/3x2x256x64x24/0T';
% resdir = '/adiabatic_electrons/5x3x192x48x24/0T';
% resdir = '/adiabatic_electrons/5x2x256x64x24/NT';
resdir = '/hot_electrons/256x64x32/0T/noise_init';
% resdir = '/hot_electrons/L_300/256x64x32/NT';
DATADIR = [PARTITION,resdir,'/'];
read_flux_out_XX(DATADIR,1,10);
%%
PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/no_gradN/';
Models = {'ion_scale/5x3x192x48x24',...
'adiabatic_electrons/5x3x192x48x24',...
'hot_electrons/256x64x32'};
% Models = {'ion_scale/5x2x256x64x24',...
% 'adiabatic_electrons/5x2x256x64x24',...
% 'hot_electrons/256x64x32'};
clrs = lines(3);
triangularities = {'NT','0T','PT'};
figure
t = tiledlayout(1,3,'TileSpacing','Compact','Padding','Compact');
for i = 1:3
nexttile
for j = 1:3
DATADIR = [PARTITION,Models{j},'/',triangularities{i},'/'];
out = read_flux_out_XX(DATADIR,0,1);
plot(out.t, out.Qxi,'color',clrs(j,:)); hold on;
end
ylim([0 200]); ylabel('$Q_{xi}$');
xlim([0 300]); xlabel('$tc_s/R$');
end
......@@ -72,6 +72,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
......@@ -104,6 +104,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
......@@ -74,6 +74,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
......@@ -72,6 +72,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
......@@ -74,6 +74,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
......@@ -109,6 +109,7 @@ 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_FVEL = 1; % Output flag for fluid velocity
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)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment