From e7ff3cfe8d39d97b58706e172dc4a520b4b021be Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 16 Apr 2024 13:05:33 +0200
Subject: [PATCH] matlab scripts update
---
matlab/compute/HermitePoly.m | 40 --------
matlab/compute/compute_fa_2D.m | 94 ++++++------------
matlab/extract_fig_data.m | 7 +-
matlab/plot/plot_metric.m | 6 --
matlab/plot/show_moments_spectrum.m | 9 +-
matlab/setup.m | 1 +
matlab/write_fort90.m | 1 +
wk/Ch7_HF_analysis.m | 17 ++--
wk/analysis_gene.m | 73 +-------------
wk/fast_analysis.m | 96 ++++++++++++-------
wk/lin_run_script.m | 41 ++++++--
wk/lin_scan_script_P_ky.m | 16 +++-
wk/load_metadata_scan.m | 36 +++++--
wk/multi_fidelity_analysis_triangularity.m | 55 +++++++++++
wk/parameters/lin_AUG.m | 1 +
wk/parameters/lin_DIIID_AUSTIN.m | 1 +
wk/parameters/lin_DIIID_LM_rho90.m | 1 +
wk/parameters/lin_DIIID_LM_rho95.m | 1 +
wk/parameters/lin_DIIID_LM_rho95_HEL.m | 1 +
wk/parameters/lin_DIIID_data.m | 1 +
wk/parameters/lin_DTT_HM_rho85.m | 1 +
wk/parameters/lin_DTT_HM_rho98.m | 1 +
wk/parameters/lin_ETG_adiab_i.m | 1 +
wk/parameters/lin_Entropy.m | 1 +
wk/parameters/lin_GASTD.m | 1 +
wk/parameters/lin_ITG.m | 1 +
wk/parameters/lin_Ivanov.m | 1 +
wk/parameters/lin_JET_rho97.m | 1 +
wk/parameters/lin_KBM.m | 1 +
wk/parameters/lin_RHT.m | 1 +
wk/parameters/lin_STEP_EC_HD_psi49.m | 1 +
wk/parameters/lin_STEP_EC_HD_psi71.m | 1 +
wk/plot_dist_func_and_trap_pass_el.m | 22 +++--
...iangularitypaper_nonlinear_conv_analysis.m | 30 ++++++
34 files changed, 304 insertions(+), 258 deletions(-)
delete mode 100644 matlab/compute/HermitePoly.m
create mode 100644 wk/multi_fidelity_analysis_triangularity.m
create mode 100644 wk/triangularitypaper_nonlinear_conv_analysis.m
diff --git a/matlab/compute/HermitePoly.m b/matlab/compute/HermitePoly.m
deleted file mode 100644
index e5c0bc58..00000000
--- a/matlab/compute/HermitePoly.m
+++ /dev/null
@@ -1,40 +0,0 @@
-
-% 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
diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index f2e9552e..4939f84f 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -1,8 +1,8 @@
-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
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index a9689411..79df88e4 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -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));
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 9a88d3b5..2d53dec0 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -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
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index f193f118..e51a2d47 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -1,6 +1,6 @@
-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';
diff --git a/matlab/setup.m b/matlab/setup.m
index ba45da52..562293dd 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -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
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 0d81ae12..01e41e08 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -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');
diff --git a/wk/Ch7_HF_analysis.m b/wk/Ch7_HF_analysis.m
index 750ec4b9..d55efe52 100644
--- a/wk/Ch7_HF_analysis.m
+++ b/wk/Ch7_HF_analysis.m
@@ -1,22 +1,21 @@
-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'};
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 15c758bd..695f5991 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -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$',...
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 4bc0c79a..490a4f92 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -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
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 79127a73..402d2caf 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -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
diff --git a/wk/lin_scan_script_P_ky.m b/wk/lin_scan_script_P_ky.m
index b0780091..2028376c 100644
--- a/wk/lin_scan_script_P_ky.m
+++ b/wk/lin_scan_script_P_ky.m
@@ -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));
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 823f506c..44fe52da 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -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,:));
diff --git a/wk/multi_fidelity_analysis_triangularity.m b/wk/multi_fidelity_analysis_triangularity.m
new file mode 100644
index 00000000..f54dfdd9
--- /dev/null
+++ b/wk/multi_fidelity_analysis_triangularity.m
@@ -0,0 +1,55 @@
+
+% 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
+
diff --git a/wk/parameters/lin_AUG.m b/wk/parameters/lin_AUG.m
index c3714ec1..35d14e86 100644
--- a/wk/parameters/lin_AUG.m
+++ b/wk/parameters/lin_AUG.m
@@ -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)
diff --git a/wk/parameters/lin_DIIID_AUSTIN.m b/wk/parameters/lin_DIIID_AUSTIN.m
index d5bc43ca..6c8492ba 100644
--- a/wk/parameters/lin_DIIID_AUSTIN.m
+++ b/wk/parameters/lin_DIIID_AUSTIN.m
@@ -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)
diff --git a/wk/parameters/lin_DIIID_LM_rho90.m b/wk/parameters/lin_DIIID_LM_rho90.m
index 85c52caf..38a65cee 100644
--- a/wk/parameters/lin_DIIID_LM_rho90.m
+++ b/wk/parameters/lin_DIIID_LM_rho90.m
@@ -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)
diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m
index ba750005..0894f58a 100644
--- a/wk/parameters/lin_DIIID_LM_rho95.m
+++ b/wk/parameters/lin_DIIID_LM_rho95.m
@@ -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)
diff --git a/wk/parameters/lin_DIIID_LM_rho95_HEL.m b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
index 0a6515a6..ff5ac727 100644
--- a/wk/parameters/lin_DIIID_LM_rho95_HEL.m
+++ b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
@@ -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)
diff --git a/wk/parameters/lin_DIIID_data.m b/wk/parameters/lin_DIIID_data.m
index 50ac9c72..a01ed51a 100644
--- a/wk/parameters/lin_DIIID_data.m
+++ b/wk/parameters/lin_DIIID_data.m
@@ -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)
diff --git a/wk/parameters/lin_DTT_HM_rho85.m b/wk/parameters/lin_DTT_HM_rho85.m
index 8fb86b3a..6e3a43ab 100644
--- a/wk/parameters/lin_DTT_HM_rho85.m
+++ b/wk/parameters/lin_DTT_HM_rho85.m
@@ -83,6 +83,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)
diff --git a/wk/parameters/lin_DTT_HM_rho98.m b/wk/parameters/lin_DTT_HM_rho98.m
index b0e204c0..59c20573 100644
--- a/wk/parameters/lin_DTT_HM_rho98.m
+++ b/wk/parameters/lin_DTT_HM_rho98.m
@@ -84,6 +84,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)
diff --git a/wk/parameters/lin_ETG_adiab_i.m b/wk/parameters/lin_ETG_adiab_i.m
index ba65082d..83920c24 100644
--- a/wk/parameters/lin_ETG_adiab_i.m
+++ b/wk/parameters/lin_ETG_adiab_i.m
@@ -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)
diff --git a/wk/parameters/lin_Entropy.m b/wk/parameters/lin_Entropy.m
index ef242a93..74021cd6 100644
--- a/wk/parameters/lin_Entropy.m
+++ b/wk/parameters/lin_Entropy.m
@@ -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)
diff --git a/wk/parameters/lin_GASTD.m b/wk/parameters/lin_GASTD.m
index 0fcb56a9..af24e3ef 100644
--- a/wk/parameters/lin_GASTD.m
+++ b/wk/parameters/lin_GASTD.m
@@ -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)
diff --git a/wk/parameters/lin_ITG.m b/wk/parameters/lin_ITG.m
index e93816e2..bd535c1d 100644
--- a/wk/parameters/lin_ITG.m
+++ b/wk/parameters/lin_ITG.m
@@ -68,6 +68,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)
diff --git a/wk/parameters/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m
index ef267364..718bc091 100644
--- a/wk/parameters/lin_Ivanov.m
+++ b/wk/parameters/lin_Ivanov.m
@@ -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)
diff --git a/wk/parameters/lin_JET_rho97.m b/wk/parameters/lin_JET_rho97.m
index 9fbd5884..28a1822f 100644
--- a/wk/parameters/lin_JET_rho97.m
+++ b/wk/parameters/lin_JET_rho97.m
@@ -105,6 +105,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)
diff --git a/wk/parameters/lin_KBM.m b/wk/parameters/lin_KBM.m
index b848aa41..fada5b9c 100644
--- a/wk/parameters/lin_KBM.m
+++ b/wk/parameters/lin_KBM.m
@@ -69,6 +69,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)
diff --git a/wk/parameters/lin_RHT.m b/wk/parameters/lin_RHT.m
index a7e2e9a3..96edd139 100644
--- a/wk/parameters/lin_RHT.m
+++ b/wk/parameters/lin_RHT.m
@@ -70,6 +70,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)
diff --git a/wk/parameters/lin_STEP_EC_HD_psi49.m b/wk/parameters/lin_STEP_EC_HD_psi49.m
index 8ca4cb51..5858b418 100644
--- a/wk/parameters/lin_STEP_EC_HD_psi49.m
+++ b/wk/parameters/lin_STEP_EC_HD_psi49.m
@@ -91,6 +91,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)
diff --git a/wk/parameters/lin_STEP_EC_HD_psi71.m b/wk/parameters/lin_STEP_EC_HD_psi71.m
index 41c60816..a99d2ce6 100644
--- a/wk/parameters/lin_STEP_EC_HD_psi71.m
+++ b/wk/parameters/lin_STEP_EC_HD_psi71.m
@@ -91,6 +91,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)
diff --git a/wk/plot_dist_func_and_trap_pass_el.m b/wk/plot_dist_func_and_trap_pass_el.m
index b07ad0bf..ba958c28 100644
--- a/wk/plot_dist_func_and_trap_pass_el.m
+++ b/wk/plot_dist_func_and_trap_pass_el.m
@@ -1,19 +1,23 @@
+[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
options.SHOW_FLUXSURF = 0;
options.SHOW_METRICS = 0;
options.SHOW_CURVOP = 0;
[fig, geo_arrays] = plot_metric(data,options);
eps_eff = 1/max((geo_arrays.hatB));
-T = [10];
-s = linspace(-3,3,64);
-x = linspace(0,8,48);
-[SS,XX,FF] = compute_fa_2D(data,'i',s,x,T);
+T = [40];
+s = linspace(-4,4,64);
+x = linspace( 0,4,64);
+[SS,XX,FF,FAM] = compute_fa_2D(data,'i',s,x,T);
+
+% FF = FF - FAM;
figure
-% contour(SS,XX,log10(FF))
-pca = pcolor(SS,XX,(FF)); set(pca,'EdgeColor','None');
-set(gca,'ColorScale','log'); %shading interp;
+contourf(SS,XX,(FF),20)
+% pca = pcolor(SS,XX,(FF)); set(pca,'EdgeColor','None');
+% set(gca,'ColorScale','log'); %shading interp;
% clim([1e-4 1]);
hold on
-plot(s,(s.^2/(1+data.fort_00.GEOMETRY.eps)),'--w');
-plot(s,eps_eff*(s.^2),'--k');
\ No newline at end of file
+% plot(s,(s.^2/(1+data.fort_00.GEOMETRY.eps)),'--w');
+% plot(s,eps_eff*(s.^2),'--k');
+colormap(bluewhitered)
\ No newline at end of file
diff --git a/wk/triangularitypaper_nonlinear_conv_analysis.m b/wk/triangularitypaper_nonlinear_conv_analysis.m
new file mode 100644
index 00000000..cf557d7b
--- /dev/null
+++ b/wk/triangularitypaper_nonlinear_conv_analysis.m
@@ -0,0 +1,30 @@
+PJ = [3*2; 5*3; 5*3; 9*5];
+Gx = [...
+ 2.8750,0.8484;...
+ 3.6025,0.8057;...
+ 3.3209,0.3283;...
+ 9.4672,12.3946;...
+ ];
+Qxe= [...
+ 42.3161,6.4356;...
+ 23.9335,5.4995;...
+ 22.9832,2.3751;...
+ 21.2910,3.1200;...
+ ];
+Qxi= [...
+ 54.7350,12.8821;...
+ 64.2177,15.8709;...
+ 60.5471,7.2580;...
+ 63.4687,10.1195;...
+ ];
+
+
+
+figure;
+errorbar(PJ,Qxe(:,1),Qxe(:,2)); hold on;
+errorbar(PJ,Qxi(:,1),Qxi(:,2));
+errorbar(PJ,Gx(:,1),Gx(:,2));
+set(gca,'Xscale','log')
+xlim([5 100])
+xticks([6 15 45]);
+xticklabels({'(2,1)','(4,2)','(8,4)'});
\ No newline at end of file
--
GitLab