diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m index 01cf310..50b3593 100644 --- a/matlab/compute/compute_fa_2D.m +++ b/matlab/compute/compute_fa_2D.m @@ -1,77 +1,107 @@ function [SS,XX,FF] = compute_fa_2D(data, options) %% 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)); % Laguerre Lj = @(j,x) polyval(LaguerrePoly(j),x); % Maxwellian factor FaM = @(s,x) exp(-s.^2-x); %meshgrid for efficient evaluation [SS, XX] = meshgrid(options.SPAR, options.XPERP); [~, ikx0] = min(abs(data.kx)); [~, iky0] = min(abs(data.ky)); +kx_ = data.kx; ky_ = data.ky; switch options.SPECIE case 'e' Napj_ = data.Nepj; parray = double(data.Pe); jarray = double(data.Je); case 'i' Napj_ = data.Nipj; parray = double(data.Pi); jarray = double(data.Ji); end switch options.Z case 'avg' Napj_ = mean(Napj_,5); + phi_ = mean(data.PHI,3); otherwise - [~,iz] = min(abs(options.Z-data.z)); + iz = options.Z; Napj_ = Napj_(:,:,:,:,iz,:); + phi_ = data.PHI(:,:,iz); end Napj_ = squeeze(Napj_); frames = options.T; for it = 1:numel(options.T) [~,frames(it)] = min(abs(options.T(it)-data.Ts5D)); end +frames = unique(frames); Napj_ = mean(Napj_(:,:,:,:,frames),5); Napj_ = squeeze(Napj_); Np = numel(parray); Nj = numel(jarray); -FF = zeros(data.Nkx,data.Nky,numel(options.XPERP),numel(options.SPAR)); -% FF = zeros(numel(options.XPERP),numel(options.SPAR)); - -FAM = FaM(SS,XX); -for ip_ = 1:Np - p_ = parray(ip_); - HH = Hp(p_,SS); +if options.non_adiab for ij_ = 1:Nj - j_ = jarray(ij_); - LL = Lj(j_,XX); -% FF = FF + Napj_(ip_,ij_,ikx0,iky0)*HH.*LL.*FAM; - HLF = HH.*LL.*FAM; for ikx = 1:data.Nkx - for iky = 1:data.Nky - FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF; + for iky = 1:data.Nky + kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2); + Napj_(1,ij_,ikx,iky) = Napj_(1,ij_,ikx,iky) + kernel(ij_,kp_)*phi_(ikx,iky); end end end end + +if options.RMS + FF = zeros(data.Nkx,data.Nky,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); + HLF = HH.*LL.*FAM; + for ikx = 1:data.Nkx + for iky = 1:data.Nky + FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF; + end + end + 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 = (FF.*conj(FF)); %|f_a|^2 % FF = abs(FF); %|f_a| -FF = sqrt(squeeze(mean(mean(FF,1),2))); %sqrt(<|f_a|>kx,ky) +if options.RMS +% FF = squeeze(mean(mean(sqrt(FF),1),2)); %sqrt(<|f_a|^2>kx,ky) + FF = sqrt(squeeze(mean(mean(FF,1),2))); %<|f_a|>kx,ky +else + FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y) +end + FF = FF./max(max(FF)); FF = FF'; -% FF = FF.*conj(FF); % FF = sqrt(FF); -% FF = FF./max(max(FF)); % FF = FF'; end \ No newline at end of file diff --git a/matlab/plot/plot_fa_gene.m b/matlab/plot/plot_fa_gene.m new file mode 100644 index 0000000..5a228a1 --- /dev/null +++ b/matlab/plot/plot_fa_gene.m @@ -0,0 +1,97 @@ +function plot_fa_gene(OPTIONS) +folder = OPTIONS.folder; +TIMES = OPTIONS.times; +specie = OPTIONS.specie; +PLT_FCT= OPTIONS.PLT_FCT; + +file = 'coord.dat.h5'; +vp = h5read([folder,file],'/coord/vp'); +mu = h5read([folder,file],'/coord/mu'); +z = h5read([folder,file],'/coord/z'); +[XX,SS] = meshgrid(mu,vp); + +switch OPTIONS.Z + case 'avg' + zcomp_name = ' z-avg'; + zcomp = @(x) squeeze(mean(x,1)); + otherwise + zcomp_name = [' z=',sprintf('%2.2f',z(OPTIONS.Z))]; + zcomp = @(x) squeeze(x(OPTIONS.Z,:,:)); +end + +[~,iv0] = min(abs(vp)); +[~,im0] = min(abs(mu)); + +file = 'vsp.dat.h5'; +time = h5read([folder,file],'/vsp/time'); + +fig = figure; + +G_t = []; +Gdata = 0; +for T = TIMES +[~, it] = min(abs(time-T)); +tmp = h5read([folder,file],['/vsp/',OPTIONS.FIELD,'/',sprintf('%10.10d',it-1)]); +Gdata = Gdata + tmp; +G_t = [G_t time(it)]; +end +Gdata = Gdata ./ numel(TIMES); + +if OPTIONS.ONED + switch specie + case 'e' + FFa = squeeze(Gdata(1,:,:,2)); + FFa = abs(FFa)./max(max(abs(FFa))); + subplot(1,2,1) + plot(vp,FFa(:,im0)); hold on; + subplot(1,2,2) + plot(mu,FFa(iv0,:)); hold on; + case 'i' + FFa = squeeze(Gdata(1,:,:,1)); + FFa = abs(FFa)./max(max(abs(FFa))); + end + + subplot(1,2,1) + plot(vp,FFa(:,im0)); hold on; + legend(specie) + xlabel('$v_\parallel, (\mu=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$'); + + subplot(1,2,2) + plot(mu,FFa(iv0,:)); hold on; + legend(specie) + xlabel('$\mu, (v_\parallel=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$'); +else +switch specie + case 'e' + name = '$f_e(v_\parallel,\mu_p)$'; + FFa = zcomp(squeeze(Gdata)); + FFa = abs(FFa)./max(max(abs(FFa))); + switch PLT_FCT + case 'contour' + contour(SS,XX,FFa,128); + case 'contourf' + pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none') + case 'pcolor' + pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp + end + case 'i' + name = '$f_i(v_\parallel,\mu_p)$'; + FFa = zcomp(squeeze(Gdata)); + FFa = abs(FFa)./max(max(abs(FFa))); + switch PLT_FCT + case 'contour' + contour(SS,XX,FFa,128); + case 'contourf' + pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none') + case 'pcolor' + pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp + end +end +if numel(TIMES) == 1 + title(['Gene ',name,zcomp_name,', $t = ',sprintf('%2.2f',time(it)),'$']); +else + title(['Gene ',name,zcomp_name,', averaged $t\in$[',sprintf('%2.2f',G_t(1)),',',sprintf('%2.2f',G_t(end)),']']); +end + +clear FFi FFe tmp Gdata +end \ No newline at end of file diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m index 6a08d4a..ed4d356 100644 --- a/wk/gene_analysis_3D.m +++ b/wk/gene_analysis_3D.m @@ -1,87 +1,181 @@ % 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.2/'; % 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/'; gene_data = load_gene_data(folder); % gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data); 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_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) options.INTERP = 1; fig = plot_radial_transport_and_spacetime(gene_data,options); % save_figure(data,fig) end if 0 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'xz'; % options.NAME ='f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [100 300 900]; gene_data.a = data.EPS * 2000; fig = photomaton(gene_data,options); save_figure(gene_data,fig) end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 0; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = 000:200; gene_data.a = data.EPS * 2000; create_film(gene_data,options,'.gif') end 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$',... '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'}; figure; subplot(311) for i = 1:6 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry'); subplot(312) for i = 7:10 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); subplot(313) for i = 11:16 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); end +% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/'; +folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/'; +% folder = '/misc/gene_results/cyclone/s_alpha_output_1.0/'; +gene_data = load_gene_data(folder); +gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data); +if 1 +%% Space time diagramm (fig 11 Ivanov 2020) +TAVG_0 = 0.8*gene_data.Ts3D(end); TAVG_1 = gene_data.Ts3D(end); % Averaging times duration +% chose your field to plot in spacetime diag (uzf,szf,Gx) +field = 'phi'; +compz = 'avg'; +nmvm = 1; +fig = plot_radial_transport_and_spacetime(gene_data,TAVG_0,TAVG_1,field,nmvm,compz); +% save_figure(data,fig) +end + +if 0 +%% 2D snapshots +% Options +options.INTERP = 1; +options.POLARPLOT = 0; +options.AXISEQUAL = 1; +options.NAME = '\phi'; +% options.NAME = 'n_i'; +% options.NAME = 'T_i'; +% options.NAME = '\Gamma_x'; +% options.NAME = 'k^2n_e'; +options.PLAN = 'xz'; +% options.NAME ='f_e'; +% options.PLAN = 'sx'; +options.COMP = 'avg'; +options.TIME = [100 300 900]; +gene_data.a = data.EPS * 2000; +fig = photomaton(gene_data,options); +save_figure(gene_data,fig) +end + +if 0 +%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Options +options.INTERP = 0; +options.POLARPLOT = 0; +options.NAME = '\phi'; +% options.NAME = 'v_y'; +% options.NAME = 'n_i^{NZ}'; +% options.NAME = '\Gamma_x'; +% options.NAME = 'n_i'; +options.PLAN = 'xy'; +% options.NAME = 'f_e'; +% options.PLAN = 'sx'; +options.COMP = 'avg'; +options.TIME = 000:170; +gene_data.a = data.EPS * 2000; +create_film(gene_data,options,'.gif') +end + +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$',... + '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'}; +figure; +subplot(311) + for i = 1:6 + plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; + end + xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry'); + +subplot(312) + for i = 7:10 + plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; + end + xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); + +subplot(313) + for i = 11:16 + plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; + end + xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); + +end + +if 0 +%% Show f_i(vpar,mu) +options.times = 200:600; +options.specie = 'i'; +options.PLT_FCT = 'pcolor'; +options.folder = folder; +options.Z = 'avg'; +options.FIELD = ''; +options.ONED = 0; +% options.FIELD = 'Q_es'; +plot_fa_gene(options); +end