diff --git a/matlab/compute/compute_fa.m b/matlab/compute/compute_fa.m
index a580a2ee399ff6db9b8102c58fe05b54d6aa2356..45bb4f9088ba93fa2b69fcd2cfaa66a0b5f38bde 100644
--- a/matlab/compute/compute_fa.m
+++ b/matlab/compute/compute_fa.m
@@ -32,7 +32,7 @@ switch options.Z
[~,iz] = min(abs(options.Z-data.z));
Napj_ = Napj_(:,:,:,:,iz,:);
end
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
frames = options.T;
for it = 1:numel(options.T)
@@ -41,7 +41,7 @@ end
Napj_ = mean(Napj_(:,:,:,:,frames),5);
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
diff --git a/matlab/compute/compute_fa_1D.m b/matlab/compute/compute_fa_1D.m
index e4774f197b5df9fea337e7953f84a6f485d4e49b..1d24967649aaf29c14351c402520a18b2996613e 100644
--- a/matlab/compute/compute_fa_1D.m
+++ b/matlab/compute/compute_fa_1D.m
@@ -38,7 +38,7 @@ switch options.iz
Napj_ = Napj_(:,:,:,:,iz,:);
phi_ = data.PHI(:,:,iz);
end
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
frames = options.T;
for it = 1:numel(options.T)
@@ -47,7 +47,7 @@ end
Napj_ = mean(Napj_(:,:,:,:,frames),5);
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
if options.non_adiab
for ij_ = 1:Nj
@@ -129,8 +129,8 @@ end
Fs = real(Fs.*conj(Fs)); %|f_a|^2
Fx = real(Fx.*conj(Fx)); %|f_a|^2
if options.RMS
-Fs = squeeze(sqrt(mean(mean(Fs,1),2))); %sqrt(<|f_a|^2>kx,ky)
-Fx = squeeze(sqrt(mean(mean(Fx,1),2))); %sqrt(<|f_a|^2>kx,ky)
+Fs = squeeze(sqrt(sum(sum(Fs,1),2))); %sqrt(<|f_a|^2>kx,ky)
+Fx = squeeze(sqrt(sum(sum(Fx,1),2))); %sqrt(<|f_a|^2>kx,ky)
end
Fs = Fs./max(max(Fs));
Fx = Fx./max(max(Fx));
diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index 29a8b093627e77d9d88722421df8d13acba1598e..a9c099318a1e3ccf7e3a17b46009e592e5beb5ba 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -35,7 +35,7 @@ switch options.iz
Napj_ = Napj_(:,:,:,:,iz,:);
phi_ = data.PHI(:,:,iz);
end
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
frames = options.T;
for it = 1:numel(options.T)
@@ -45,7 +45,7 @@ frames = unique(frames);
Napj_ = mean(Napj_(:,:,:,:,frames),5);
-Napj_ = squeeze(Napj_);
+% Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index 278fe024808b5f8264211c536e4c964bf3d1f4c3..fddaa5e60ca7a56b157e4f0d7157814494ac26e4 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -17,11 +17,7 @@ TNAME = [];
FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
for i_ = 1:numel(FRAMES)
subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
- if ~strcmp(OPTIONS.PLAN,'kxky')
- scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize
- else
- scale = 1;
- end
+ scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize
if ~strcmp(OPTIONS.PLAN,'sx')
tshot = DATA.Ts3D(FRAMES(i_));
pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 7e31997c422d39c062afb8f6037a127e6c2fd809..968f87656a430a605a3f37158b90c7c0b0c747e5 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -10,8 +10,11 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
+ disp(['G_x=',sprintf('%2.2f',Gx_infty_avg),'+-',sprintf('%2.2f',Gx_infty_std)]);
+ disp(['Q_x=',sprintf('%2.2f',Qx_infty_avg),'+-',sprintf('%2.2f',Qx_infty_std)]);
f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
+ ikzf = min([ikzf,DATA.Nky]);
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
@@ -51,20 +54,22 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
subplot(311)
- yyaxis left
- plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
- plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
- plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
- 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
- ylabel('$\Gamma_x$')
- ylim([0,5*abs(Gx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
- yyaxis right
+% yyaxis left
+% plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+% plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
+% plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
+% 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
+% ylabel('$\Gamma_x$')
+% ylim([0,5*abs(Gx_infty_avg)]);
+% xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+% yyaxis right
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
- plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
+% plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
ylabel('$Q_x$')
- ylim([0,5*abs(Qx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+ ylim([0,2*abs(Qx_infty_avg)]);
+ xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index 67ed61da2d0e3ff16c8ca3eb78b990951c759788..eb284c69a1d2333e5a792135a60e932b57f46ae1 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -1,9 +1,11 @@
+
function [ fig ] = statistical_transport_averaging( data, options )
-scale = (1/data.Nx/data.Ny)^2;
+scale = data.scale;
Trange = options.T;
[~,it0] = min(abs(Trange(1)-data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D));
gamma = data.PGAMMA_RI(it0:it1)*scale;
+Qx = data.HFLUX_X(it0:it1)*scale;
dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
% if ~dt_const
% disp('DT not const on given interval');
@@ -22,13 +24,22 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0));
fig = figure;
-% subplot(211)
+ subplot(211)
+ plot(time_seg,transp_seg_avg,'-'); hold on;
+ xlabel('Averaging time'); ylabel('$\langle\Gamma_x\rangle_{\tau}$');
+ legend(['$\Gamma_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))])
+ title(sprintf('Transport averaging from t=%2.2f',data.Ts0D(it0)));
+
+ for Np = 1:Ntot % Loop on the number of segments
+ transp_seg_avg(Np) = mean(Qx(1:Np));
+ end
+
+
+ subplot(212)
+
plot(time_seg,transp_seg_avg,'-'); hold on;
- xlabel('Averaging time'); ylabel('transport value');
+ xlabel('Averaging time'); ylabel('$\langle Q_x\rangle_{\tau}$');
+ legend(['$Q_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))])
-% subplot(212)
-% errorbar(N_seg,transp_seg_avg,transp_seg_std);
-% xlabel('Averaging #points'); ylabel('transport value');
-% end
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 201d0a49725f55652ce61b43896c42e71daa6625..b6409f70188b3385a1b169b7d5238b90cf8e0ec6 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -447,6 +447,7 @@ switch OPTIONS.NAME
for it = 1:numel(OPTIONS.TIME)
[~,it0_] =min(abs(OPTIONS.TIME(it)-DATA.Ts5D));
OPTIONS.T = DATA.Ts5D(it0_);
+ OPTIONS.Z = OPTIONS.COMP;
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
case 'f_e'
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index fbe4b7795798620e9833880d5a4e189453673ef5..97d8414a775420cf983d4006ac30142eee623687 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -22,7 +22,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.98*data.Ts3D(end);
+options.TAVG_0 = 0.5*data.Ts3D(end); data.scale = (1/data.Nx/data.Ny)^2;
options.TAVG_1 = data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -34,14 +34,14 @@ end
if 0
%% statistical transport averaging
-options.T = [16000 17000];
+options.T = [200 500];
fig = statistical_transport_averaging(data,options);
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
@@ -50,11 +50,12 @@ options.NAME = '\phi';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xz';
-% options.NAME = 'f_e';
+% options.NAME = 'f_i';
% options.PLAN = 'sx';
-options.COMP = 9;
-% options.TIME = dat.Ts5D;
-options.TIME = [0:1:2000];
+options.COMP = 'avg';
+% options.TIME = data.Ts5D(end-30:end);
+options.TIME = data.Ts3D;
+% options.TIME = [700:1100];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -72,14 +73,14 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
-% options.NAME = 'f_i';
-% options.PLAN = 'sx';
+% options.PLAN = 'kxky';
+options.NAME = 'f_i';
+options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [1600 1800 2000];
+options.TIME = [200 400 700];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
-save_figure(data,fig)
+% save_figure(data,fig)
end
if 0
@@ -100,26 +101,26 @@ if 0
% options.XPERP = linspace( 0,6,64);
options.SPAR = gene_data.vp';
options.XPERP = gene_data.mu';
-options.iz = 9;
-options.T = 30;
+options.iz = 'avg';
+options.T = [500];
options.PLT_FCT = 'contour';
options.ONED = 0;
-options.non_adiab = 1;
+options.non_adiab = 0;
options.SPECIE = 'i';
options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
fig = plot_fa(data,options);
-save_figure(data,fig)
+% save_figure(data,fig)
end
if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
-options.P2J = 1;
-options.ST = 0;
+options.P2J = 0;
+options.ST = 1;
options.PLOT_TYPE = 'space-time';
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
options.JOBNUM = 0;
-options.TIME = [50];
+options.TIME = [300:500];
options.specie = 'i';
options.compz = 'avg';
fig = show_moments_spectrum(data,options);
@@ -129,7 +130,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = 1000:1200;
+options.TIME = 650:800;
options.NORM =1;
options.NAME = '\phi';
% options.NAME = 'n_i';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 86e064c561aa6b39a9b9086478f451dd12a62ed5..c112903b01a5ee6c783fee1908efd8a5708a5d0d 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -1,17 +1,19 @@
% 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/s_alpha_output_0.8/';
+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/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
+% folder = '/misc/gene_results/LD_zpinch_1.6/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.8*gene_data.Ts3D(end);
+options.TAVG_0 = 0.2*gene_data.Ts3D(end);
options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
@@ -21,6 +23,12 @@ fig = plot_radial_transport_and_spacetime(gene_data,options);
save_figure(gene_data,fig)
end
+if 0
+%% statistical transport averaging
+options.T = [100 500];
+fig = statistical_transport_averaging(gene_data,options);
+end
+
if 0
%% 2D snapshots
% Options
@@ -32,11 +40,11 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'xz';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [500];
+options.TIME = [0];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig)
@@ -52,7 +60,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'xz';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -95,11 +103,28 @@ options.PLT_FCT = 'contour';
options.folder = folder;
options.iz = 9;
options.FIELD = '<f_>';
-options.ONED = 0;
+options.ONED = 1;
% options.FIELD = 'Q_es';
plot_fa_gene(options);
end
+if 0
+%% Time averaged spectrum
+options.TIME = 300:600;
+options.NORM =1;
+options.NAME = '\phi';
+% options.NAME = 'n_i';
+% options.NAME ='\Gamma_x';
+options.PLAN = 'kxky';
+options.COMPZ = 'avg';
+options.OK = 0;
+options.COMPXY = 'avg';
+options.COMPT = 'avg';
+options.PLOT = 'semilogy';
+fig = spectrum_1D(data,options);
+% save_figure(data,fig)
+end
+
if 0
%% Mode evolution
options.NORMALIZED = 0;
diff --git a/wk/cbc_qx_results.m b/wk/cbc_qx_results.m
new file mode 100644
index 0000000000000000000000000000000000000000..342175f9e5123efec56aad95b4812f2d435da25b
--- /dev/null
+++ b/wk/cbc_qx_results.m
@@ -0,0 +1,18 @@
+cbc = [0080 0100 0120];
+gm42 = [15.1 30.2 51.2];
+gm42_err = [2.00 04.8 05.3];
+gm84 = [6.74 25.6 39.4];
+gm84_err = [0.00 00.0 00.0];
+gne = [8.44 24.1 43.5];
+gne_err = [1.55 04.3 05.9]/2;
+kN = [1.776 2.22 2.664];
+kT = [5.568 6.96 8.352];
+figure
+errorbar(cbc,gm42,gm42_err,'o-','LineWidth',1.5); hold on;
+errorbar(cbc,gm84,gm84_err,'o-','LineWidth',1.5);
+errorbar(cbc,gne,gne_err,'x-k','LineWidth',1.5);
+
+set(gca, 'YScale', 'log')
+
+legend('GM (4,2)','GM (8,4)','Gene')
+xlabel('CBC drive [\%]'); ylabel('Radial Heat Flux $Q_x^\infty$');
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index b028f6c0c52db2bcff1b088fefeb07747b8f3b52..0b3663bf923b35091f0a5cbec5e964747547d644 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -4,17 +4,12 @@ helazdir = '/home/ahoffman/HeLaZ/';
% if 1% Local results
outfile ='';
outfile ='';
-outfile ='';
-% outfile ='shearless_cyclone/128x128x16x5x3_start';
-% outfile ='shearless_cyclone/128x128x16_CBC_100';
-outfile ='shearless_cyclone/128x128x16_CBC_120';
-% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
-% outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
-% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
-% outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
-% outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
-% outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
-% outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/';
-% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/';
-% outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1';
+% outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
+% outfile ='shearless_cyclone/64x32x16x5x3_CBC_100';
+% outfile ='shearless_cyclone/64x32x16x5x3_CBC_120';
+
+% outfile ='shearless_cyclone/64x32x16x9x5_CBC_080';
+% outfile ='shearless_cyclone/64x32x16x9x5_CBC_100';
+% outfile ='shearless_cyclone/64x32x16x9x5_CBC_120';
+
run analysis_HeLaZ
diff --git a/wk/heat_flux_postproc.m b/wk/heat_flux_postproc.m
new file mode 100644
index 0000000000000000000000000000000000000000..500ebf9a085650f3453e11fe8df57dccdd03c92e
--- /dev/null
+++ b/wk/heat_flux_postproc.m
@@ -0,0 +1,105 @@
+% Hflux_x = 0;
+% Hflux_x = 0 * data.Ts5D;
+filename = '/home/ahoffman/HeLaZ/results/shearless_cyclone/64x32x16x5x3_CBC_100/outputs_00.h5';
+kernel_i = h5read(filename,'/data/metric/kernel_i');
+Jacobian = h5read(filename,'/data/metric/Jacobian');
+STEPS = 1:numel(data.Ts5D);
+Hflux_x = 1:numel(STEPS);
+its_ = 1;
+for it = STEPS
+ t = data.Ts5D(it);
+ [~,it5d] = min(abs(data.Ts5D-t));
+ [~,it3d] = min(abs(data.Ts3D-t));
+
+ Nj = data.Jmaxi;
+ Nz = data.Nz; z = data.z; dz = data.z(2)-data.z(1);
+ kx = data.kx; ky = data.ky; Nkx = data.Nkx; Nky = data.Nky;
+ Jz = squeeze(Jacobian(:,1));
+ invjac = 1/2/pi/data.Q0;
+
+ % Factors for sum kern mom
+ c2n = @(n_) 0.5*sqrt(2);
+ c0n = @(n_) 2*n_+ 1.5;
+% c0n = @(n_) 2*n_+ 2/3;
+ c0np1 = @(n_) -(n_+1);
+ c0nm1 = @(n_) -n_;
+
+ % Factors for correction operator
+% dn = @(n_) -2*(n_+ 1.5);
+ dn = @(n_) -2*n_+ 1.5;
+ dnp1 = @(n_) (n_+1);
+ dnm1 = @(n_) n_;
+
+ % BUILD TERM TO SUM
+ sumkernmom = zeros(Nky,Nkx,Nz);
+ correct_op = zeros(Nky,Nkx,Nz);
+ for in = 1:Nj
+ n = in-1;
+ Kn = squeeze(kernel_i(in,:,:,:,1));
+ N2n = squeeze(data.Nipj(3,in ,:,:,:,it5d));
+ N0n = squeeze(data.Nipj(1,in ,:,:,:,it5d));
+ sumkernmom = sumkernmom + ...
+ Kn.* (c2n(n) .* N2n + c0n(n) .* N0n);
+
+ correct_op = correct_op + ...
+ Kn.* dn(n) .* Kn;
+
+ if(in > 1)
+ N0nm1 = squeeze(data.Nipj(1,in-1,:,:,:,it5d));
+ sumkernmom = sumkernmom + Kn.* c0nm1(n).*N0nm1;
+
+ Knm1 = squeeze(kernel_i(in-1,:,:,:,1));
+ correct_op = correct_op + Kn.* dnm1(n) .* Knm1;
+ end
+ if(in<Nj)
+ N0np1 = squeeze(data.Nipj(1,in+1,:,:,:,it5d));
+ sumkernmom = sumkernmom + Kn.* c0np1(n).*N0np1;
+
+ Knp1 = squeeze(kernel_i(in+1,:,:,:,1));
+ correct_op = correct_op + Kn.* dnp1(n) .* Knp1;
+ end
+ end
+ [~,KYY,ZZZ] = meshgrid(kx,ky,z);
+ % -- adding correction term
+ correction_term = data.PHI(:,:,:,it3d) .* correct_op/ data.scale;
+
+ % -- summing up
+ u = -1i*KYY.*data.PHI(:,:,:,it3d);
+ n = sumkernmom + correction_term;
+ clear sumkernmom correction_term correct_op KYY;
+
+ sum_kxky = conj(u).*n;
+ half_plane= sum_kxky;
+ Ny = 2*Nky-1;
+ sum_kxky = zeros(Ny,Nkx,Nz);
+ sum_kxky(1:Nky,:,:)=half_plane(:,:,:);
+ sum_kxky((Nky+1):(Ny),1,:)=conj(half_plane(Nky:-1:2,1,:));
+ sum_kxky((Nky+1):(Ny),2:Nkx,:)=conj(half_plane(Nky:-1:2,Nkx:-1:2,:));
+
+ sum_kxky = squeeze(sum(sum(sum_kxky,1),2));
+ % Z average
+ J_= 0;
+ q_ = 0;
+ for iz = 1:Nz
+% add = u(1,1,iz)*n(1,1,iz)...
+% +2*sum(real(u(2:Nky,1,iz).*conj(n(2:Nky,1,iz))),1)...
+% +sum(2*real(u(1,2:Nkx/2+1,iz).*conj(n(1,2:Nkx/2+1,iz))) ...
+% +sum(2*real(u(2:Nky,2:Nkx/2+1,iz).*conj(n(2:Nky,2:Nkx/2+1,iz)))+...
+% 2*real(u(2:Nky,Nkx/2+1:Nkx,iz).*conj(n(2:Nky,Nkx/2+1:Nkx,iz))),1),2);
+ add = sum_kxky(iz);
+ if(mod(iz,2)==1)
+ q_ = q_ + 4 * Jz(iz)*add;
+ J_ = J_ + 4 * Jz(iz);
+ else
+ q_ = q_ + 2 * Jz(iz)*add;
+ J_ = J_ + 2 * Jz(iz);
+ end
+ end
+ q_ = q_/J_;
+ Hflux_x(its_) = real(q_);
+ its_ = its_ + 1;
+end
+%%
+figure
+plot(data.Ts0D,data.HFLUX_X.*data.scale); hold on;
+plot(data.Ts5D(STEPS),Hflux_x.*data.scale)
\ No newline at end of file
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 77df8c3ca7ea9d25761489d8ea8812d317e92857..5068eff7960f325dbbf7f28e2cd35bf4ef3fffe6 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -13,7 +13,7 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.0; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 0;%2.22; % Density gradient drive
K_T = 0;%6.96; % Temperature '''
@@ -21,10 +21,10 @@ K_E = 0.0; % Electrostat '''
SIGMA_E = 0.05196152422706632;%0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
%% GRID PARAMETERS
-PMAXE = 6; % Hermite basis size of electrons
-JMAXE = 3; % Laguerre "
-PMAXI = 6; % " ions
-JMAXI = 3; % "
+PMAXE = 8; % Hermite basis size of electrons
+JMAXE = 4; % Laguerre "
+PMAXI = 8; % " ions
+JMAXI = 4; % "
NX = 2; % real space x-gridpoints
NY = 1; % '' y-gridpoints
LX = 100; % Size of the squared frequency domain
@@ -39,8 +39,8 @@ Q0 = 1.4; % safety factor
SHEAR = 0.0; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 30; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 50; % Maximal time unit
+DT = 2e-2; % Time step
SPS0D = 20; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 20; % Sampling per time unit for 2D arrays
@@ -165,10 +165,11 @@ end
if 1
%% RH TEST
-ikx = 2;
-plt = @(x) squeeze(mean(real(x(1,ikx,:,:)),3))./squeeze(mean(real(x(1,ikx,:,1)),3));
+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, plt(data.PHI));
+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