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

matlab scripts update

parent 46c46efc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 238 additions and 80 deletions
matlab/compute/compute_fa.m
+ 2
− 2
View file @ 0f50fdc7
...@@ -32,7 +32,7 @@ switch options.Z ...@@ -32,7 +32,7 @@ switch options.Z
[~,iz] = min(abs(options.Z-data.z)); [~,iz] = min(abs(options.Z-data.z));
Napj_ = Napj_(:,:,:,:,iz,:); Napj_ = Napj_(:,:,:,:,iz,:);
end end
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
frames = options.T; frames = options.T;
for it = 1:numel(options.T) for it = 1:numel(options.T)
...@@ -41,7 +41,7 @@ end ...@@ -41,7 +41,7 @@ end
Napj_ = mean(Napj_(:,:,:,:,frames),5); Napj_ = mean(Napj_(:,:,:,:,frames),5);
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray); Np = numel(parray); Nj = numel(jarray);
......
matlab/compute/compute_fa_1D.m
+ 4
− 4
View file @ 0f50fdc7
...@@ -38,7 +38,7 @@ switch options.iz ...@@ -38,7 +38,7 @@ switch options.iz
Napj_ = Napj_(:,:,:,:,iz,:); Napj_ = Napj_(:,:,:,:,iz,:);
phi_ = data.PHI(:,:,iz); phi_ = data.PHI(:,:,iz);
end end
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
frames = options.T; frames = options.T;
for it = 1:numel(options.T) for it = 1:numel(options.T)
...@@ -47,7 +47,7 @@ end ...@@ -47,7 +47,7 @@ end
Napj_ = mean(Napj_(:,:,:,:,frames),5); Napj_ = mean(Napj_(:,:,:,:,frames),5);
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
if options.non_adiab if options.non_adiab
for ij_ = 1:Nj for ij_ = 1:Nj
...@@ -129,8 +129,8 @@ end ...@@ -129,8 +129,8 @@ end
Fs = real(Fs.*conj(Fs)); %|f_a|^2 Fs = real(Fs.*conj(Fs)); %|f_a|^2
Fx = real(Fx.*conj(Fx)); %|f_a|^2 Fx = real(Fx.*conj(Fx)); %|f_a|^2
if options.RMS if options.RMS
Fs = squeeze(sqrt(mean(mean(Fs,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(mean(mean(Fx,1),2))); %sqrt(<|f_a|^2>kx,ky) Fx = squeeze(sqrt(sum(sum(Fx,1),2))); %sqrt(<|f_a|^2>kx,ky)
end end
Fs = Fs./max(max(Fs)); Fs = Fs./max(max(Fs));
Fx = Fx./max(max(Fx)); Fx = Fx./max(max(Fx));
......
matlab/compute/compute_fa_2D.m
+ 2
− 2
View file @ 0f50fdc7
...@@ -35,7 +35,7 @@ switch options.iz ...@@ -35,7 +35,7 @@ switch options.iz
Napj_ = Napj_(:,:,:,:,iz,:); Napj_ = Napj_(:,:,:,:,iz,:);
phi_ = data.PHI(:,:,iz); phi_ = data.PHI(:,:,iz);
end end
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
frames = options.T; frames = options.T;
for it = 1:numel(options.T) for it = 1:numel(options.T)
...@@ -45,7 +45,7 @@ frames = unique(frames); ...@@ -45,7 +45,7 @@ frames = unique(frames);
Napj_ = mean(Napj_(:,:,:,:,frames),5); Napj_ = mean(Napj_(:,:,:,:,frames),5);
Napj_ = squeeze(Napj_); % Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray); Np = numel(parray); Nj = numel(jarray);
......
matlab/plot/photomaton.m
+ 1
− 5
View file @ 0f50fdc7
...@@ -17,11 +17,7 @@ TNAME = []; ...@@ -17,11 +17,7 @@ TNAME = [];
FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows]) FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
for i_ = 1:numel(FRAMES) for i_ = 1:numel(FRAMES)
subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))]; 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
scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize
else
scale = 1;
end
if ~strcmp(OPTIONS.PLAN,'sx') if ~strcmp(OPTIONS.PLAN,'sx')
tshot = DATA.Ts3D(FRAMES(i_)); tshot = DATA.Ts3D(FRAMES(i_));
pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none'); pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
......
matlab/plot/plot_radial_transport_and_spacetime.m
+ 15
− 10
View file @ 0f50fdc7
...@@ -10,8 +10,11 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS) ...@@ -10,8 +10,11 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE; Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE; Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (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)); f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3))); [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
ikzf = min([ikzf,DATA.Nky]);
Ns3D = numel(DATA.Ts3D); Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky); [KX, KY] = meshgrid(DATA.kx, DATA.ky);
...@@ -51,20 +54,22 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS) ...@@ -51,20 +54,22 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
mvm = @(x) movmean(x,OPTIONS.NMVA); mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600]) FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
subplot(311) subplot(311)
yyaxis left % 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.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(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',... % 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)]); % 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
ylabel('$\Gamma_x$') % ylabel('$\Gamma_x$')
ylim([0,5*abs(Gx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]); % ylim([0,5*abs(Gx_infty_avg)]);
yyaxis right % 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.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',... 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)]); 'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
ylabel('$Q_x$') 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',[]); grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,... title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]}); ['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
......
matlab/plot/statistical_transport_averaging.m
+ 18
− 7
View file @ 0f50fdc7
function [ fig ] = statistical_transport_averaging( data, options ) function [ fig ] = statistical_transport_averaging( data, options )
scale = (1/data.Nx/data.Ny)^2; scale = data.scale;
Trange = options.T; Trange = options.T;
[~,it0] = min(abs(Trange(1)-data.Ts0D)); [~,it0] = min(abs(Trange(1)-data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D)); [~,it1] = min(abs(Trange(end)-data.Ts0D));
gamma = data.PGAMMA_RI(it0:it1)*scale; 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; dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
% if ~dt_const % if ~dt_const
% disp('DT not const on given interval'); % disp('DT not const on given interval');
...@@ -22,13 +24,22 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1; ...@@ -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)); time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0));
fig = figure; 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; 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 end
matlab/process_field.m
+ 1
− 0
View file @ 0f50fdc7
...@@ -447,6 +447,7 @@ switch OPTIONS.NAME ...@@ -447,6 +447,7 @@ switch OPTIONS.NAME
for it = 1:numel(OPTIONS.TIME) for it = 1:numel(OPTIONS.TIME)
[~,it0_] =min(abs(OPTIONS.TIME(it)-DATA.Ts5D)); [~,it0_] =min(abs(OPTIONS.TIME(it)-DATA.Ts5D));
OPTIONS.T = DATA.Ts5D(it0_); OPTIONS.T = DATA.Ts5D(it0_);
OPTIONS.Z = OPTIONS.COMP;
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS); [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end end
case 'f_e' case 'f_e'
......
wk/analysis_HeLaZ.m
+ 22
− 21
View file @ 0f50fdc7
...@@ -22,7 +22,7 @@ FMT = '.fig'; ...@@ -22,7 +22,7 @@ FMT = '.fig';
if 1 if 1
%% Space time diagramm (fig 11 Ivanov 2020) %% 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.TAVG_1 = data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces 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) % 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 ...@@ -34,14 +34,14 @@ end
if 0 if 0
%% statistical transport averaging %% statistical transport averaging
options.T = [16000 17000]; options.T = [200 500];
fig = statistical_transport_averaging(data,options); fig = statistical_transport_averaging(data,options);
end end
if 0 if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options % Options
options.INTERP = 1; options.INTERP = 0;
options.POLARPLOT = 0; options.POLARPLOT = 0;
options.NAME = '\phi'; options.NAME = '\phi';
% options.NAME = 'N_i^{00}'; % options.NAME = 'N_i^{00}';
...@@ -50,11 +50,12 @@ options.NAME = '\phi'; ...@@ -50,11 +50,12 @@ options.NAME = '\phi';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
options.PLAN = 'xz'; options.PLAN = 'xz';
% options.NAME = 'f_e'; % options.NAME = 'f_i';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 9; options.COMP = 'avg';
% options.TIME = dat.Ts5D; % options.TIME = data.Ts5D(end-30:end);
options.TIME = [0:1:2000]; options.TIME = data.Ts3D;
% options.TIME = [700:1100];
data.EPS = 0.1; data.EPS = 0.1;
data.a = data.EPS * 2000; data.a = data.EPS * 2000;
create_film(data,options,'.gif') create_film(data,options,'.gif')
...@@ -72,14 +73,14 @@ options.NAME = '\phi'; ...@@ -72,14 +73,14 @@ options.NAME = '\phi';
% options.NAME = 'T_i'; % options.NAME = 'T_i';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e'; % options.NAME = 'k^2n_e';
options.PLAN = 'xy'; % options.PLAN = 'kxky';
% options.NAME = 'f_i'; options.NAME = 'f_i';
% options.PLAN = 'sx'; options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
options.TIME = [1600 1800 2000]; options.TIME = [200 400 700];
data.a = data.EPS * 2e3; data.a = data.EPS * 2e3;
fig = photomaton(data,options); fig = photomaton(data,options);
save_figure(data,fig) % save_figure(data,fig)
end end
if 0 if 0
...@@ -100,26 +101,26 @@ if 0 ...@@ -100,26 +101,26 @@ if 0
% options.XPERP = linspace( 0,6,64); % options.XPERP = linspace( 0,6,64);
options.SPAR = gene_data.vp'; options.SPAR = gene_data.vp';
options.XPERP = gene_data.mu'; options.XPERP = gene_data.mu';
options.iz = 9; options.iz = 'avg';
options.T = 30; options.T = [500];
options.PLT_FCT = 'contour'; options.PLT_FCT = 'contour';
options.ONED = 0; options.ONED = 0;
options.non_adiab = 1; options.non_adiab = 0;
options.SPECIE = 'i'; options.SPECIE = 'i';
options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
fig = plot_fa(data,options); fig = plot_fa(data,options);
save_figure(data,fig) % save_figure(data,fig)
end end
if 0 if 0
%% Hermite-Laguerre spectrum %% Hermite-Laguerre spectrum
% options.TIME = 'avg'; % options.TIME = 'avg';
options.P2J = 1; options.P2J = 0;
options.ST = 0; options.ST = 1;
options.PLOT_TYPE = 'space-time'; options.PLOT_TYPE = 'space-time';
options.NORMALIZED = 1; options.NORMALIZED = 0;
options.JOBNUM = 0; options.JOBNUM = 0;
options.TIME = [50]; options.TIME = [300:500];
options.specie = 'i'; options.specie = 'i';
options.compz = 'avg'; options.compz = 'avg';
fig = show_moments_spectrum(data,options); fig = show_moments_spectrum(data,options);
...@@ -129,7 +130,7 @@ end ...@@ -129,7 +130,7 @@ end
if 0 if 0
%% Time averaged spectrum %% Time averaged spectrum
options.TIME = 1000:1200; options.TIME = 650:800;
options.NORM =1; options.NORM =1;
options.NAME = '\phi'; options.NAME = '\phi';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
......
wk/analysis_gene.m
+ 31
− 6
View file @ 0f50fdc7
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/'; % 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/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/linear_s_alpha_CBC_100/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/'; % 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_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/'; % 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_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_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 = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data); gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1 if 1
%% Space time diagramm (fig 11 Ivanov 2020) %% 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.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces 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) 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); ...@@ -21,6 +23,12 @@ fig = plot_radial_transport_and_spacetime(gene_data,options);
save_figure(gene_data,fig) save_figure(gene_data,fig)
end end
if 0
%% statistical transport averaging
options.T = [100 500];
fig = statistical_transport_averaging(gene_data,options);
end
if 0 if 0
%% 2D snapshots %% 2D snapshots
% Options % Options
...@@ -32,11 +40,11 @@ options.NAME = '\phi'; ...@@ -32,11 +40,11 @@ options.NAME = '\phi';
% options.NAME = 'T_i'; % options.NAME = 'T_i';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e'; % options.NAME = 'k^2n_e';
options.PLAN = 'xy'; options.PLAN = 'xz';
% options.NAME ='f_e'; % options.NAME ='f_e';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
options.TIME = [500]; options.TIME = [0];
gene_data.a = data.EPS * 2000; gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options); fig = photomaton(gene_data,options);
save_figure(gene_data,fig) save_figure(gene_data,fig)
...@@ -52,7 +60,7 @@ options.NAME = '\phi'; ...@@ -52,7 +60,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}'; % options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
options.PLAN = 'xy'; options.PLAN = 'xz';
% options.NAME = 'f_e'; % options.NAME = 'f_e';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
...@@ -95,11 +103,28 @@ options.PLT_FCT = 'contour'; ...@@ -95,11 +103,28 @@ options.PLT_FCT = 'contour';
options.folder = folder; options.folder = folder;
options.iz = 9; options.iz = 9;
options.FIELD = '<f_>'; options.FIELD = '<f_>';
options.ONED = 0; options.ONED = 1;
% options.FIELD = 'Q_es'; % options.FIELD = 'Q_es';
plot_fa_gene(options); plot_fa_gene(options);
end 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 if 0
%% Mode evolution %% Mode evolution
options.NORMALIZED = 0; options.NORMALIZED = 0;
......
wk/cbc_qx_results.m 0 → 100644
+ 18
− 0
View file @ 0f50fdc7
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
wk/header_3D_results.m
+ 8
− 13
View file @ 0f50fdc7
...@@ -4,17 +4,12 @@ helazdir = '/home/ahoffman/HeLaZ/'; ...@@ -4,17 +4,12 @@ helazdir = '/home/ahoffman/HeLaZ/';
% if 1% Local results % if 1% Local results
outfile =''; outfile ='';
outfile =''; outfile ='';
outfile =''; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
% outfile ='shearless_cyclone/128x128x16x5x3_start'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_100';
% outfile ='shearless_cyclone/128x128x16_CBC_100'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_120';
outfile ='shearless_cyclone/128x128x16_CBC_120';
% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_080';
% outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_100';
% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_120';
% 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';
run analysis_HeLaZ run analysis_HeLaZ
wk/heat_flux_postproc.m 0 → 100644
+ 105
− 0
View file @ 0f50fdc7
% 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
wk/quick_run.m
+ 11
− 10
View file @ 0f50fdc7
...@@ -13,7 +13,7 @@ EXECNAME = 'helaz3'; ...@@ -13,7 +13,7 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS %% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency NU = 0.0; % Collision frequency
TAU = 1.0; % e/i temperature ratio TAU = 1.0; % e/i temperature ratio
K_N = 0;%2.22; % Density gradient drive K_N = 0;%2.22; % Density gradient drive
K_T = 0;%6.96; % Temperature ''' K_T = 0;%6.96; % Temperature '''
...@@ -21,10 +21,10 @@ K_E = 0.0; % Electrostat ''' ...@@ -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) 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 KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
%% GRID PARAMETERS %% GRID PARAMETERS
PMAXE = 6; % Hermite basis size of electrons PMAXE = 8; % Hermite basis size of electrons
JMAXE = 3; % Laguerre " JMAXE = 4; % Laguerre "
PMAXI = 6; % " ions PMAXI = 8; % " ions
JMAXI = 3; % " JMAXI = 4; % "
NX = 2; % real space x-gridpoints NX = 2; % real space x-gridpoints
NY = 1; % '' y-gridpoints NY = 1; % '' y-gridpoints
LX = 100; % Size of the squared frequency domain LX = 100; % Size of the squared frequency domain
...@@ -39,8 +39,8 @@ Q0 = 1.4; % safety factor ...@@ -39,8 +39,8 @@ Q0 = 1.4; % safety factor
SHEAR = 0.0; % magnetic shear (Not implemented yet) SHEAR = 0.0; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS %% TIME PARMETERS
TMAX = 30; % Maximal time unit TMAX = 50; % Maximal time unit
DT = 1e-2; % Time step DT = 2e-2; % Time step
SPS0D = 20; % Sampling per time unit for 2D arrays SPS0D = 20; % Sampling per time unit for 2D arrays
SPS2D = 0; % 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 SPS3D = 20; % Sampling per time unit for 2D arrays
...@@ -165,10 +165,11 @@ end ...@@ -165,10 +165,11 @@ end
if 1 if 1
%% RH TEST %% RH TEST
ikx = 2; ikx = 2; t0 = 0; t1 = data.Ts3D(end);
plt = @(x) squeeze(mean(real(x(1,ikx,:,:)),3))./squeeze(mean(real(x(1,ikx,:,1)),3)); [~, 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 figure
plot(data.Ts3D, plt(data.PHI)); plot(data.Ts3D(it0:it1), plt(data.PHI));
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx))) title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
end end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment