Skip to content
Snippets Groups Projects
Commit 968d0648 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann 🌱
Browse files

update scripts

parent 86be0681
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/compute/compute_fa_2D_spar_mu.m 0 → 100644
+ 66
− 0
View file @ 968d0648
function [Spar,Mu,FF,FAM] = compute_fa_2D_spar_mu(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) polyval(HermitePhys_norm(p),s);
% Laguerre
Lj = @(j,x) polyval(LaguerrePoly(j),x);
% Maxwellian factor
FaM = @(s,x) exp(-s.^2-x);
%meshgrid for efficient evaluation
[Spar, Mu] = meshgrid(s, x);
switch species
case 'e'
Napj_ = data.Napjz(2,:,:,:,:);
case 'i'
Napj_ = data.Napjz(1,:,:,:,:);
end
parray = double(data.grids.Parray);
jarray = double(data.grids.Jarray);
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)
[~,frames(it)] = min(abs(T(it)-data.Ts3D));
end
frames = unique(frames);
Napj_ = mean(Napj_(:,:,:,:,frames),5);
Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
FF = zeros(numel(x),numel(s));
FAM = FaM(Spar,Mu);
for ip_ = 1:Np
p_ = parray(ip_);
HH = Hp(p_,Spar);
for ij_ = 1:Nj
j_ = jarray(ij_);
LL = Lj(j_,Mu);
HLF = HH.*LL.*FAM;
FF = FF + Napj_(ip_,ij_)*HLF;
end
end
FF = (FF.*conj(FF)); %|f_a|^2
FF = sqrt(FF); %<|f_a|>kx,ky
FF = FF./max(FF(:));
end
\ No newline at end of file
matlab/compute/compute_fa_2D_spar_sper.m 0 → 100644
+ 66
− 0
View file @ 968d0648
function [Spar,Sper,FF,FAM] = compute_fa_2D_spar_sper(data, species, spar, sper, 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) polyval(HermitePhys_norm(p),s);
% Laguerre
Lj = @(j,x) polyval(LaguerrePoly(j),x);
% Maxwellian factor
FaM = @(s,x) exp(-s.^2-x);
%meshgrid for efficient evaluation
[Spar, Sper] = meshgrid(spar, sper);
switch species
case 'e'
Napj_ = data.Napjz(2,:,:,:,:);
case 'i'
Napj_ = data.Napjz(1,:,:,:,:);
end
parray = double(data.grids.Parray);
jarray = double(data.grids.Jarray);
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)
[~,frames(it)] = min(abs(T(it)-data.Ts3D));
end
frames = unique(frames);
Napj_ = mean(Napj_(:,:,:,:,frames),5);
Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
FF = zeros(numel(sper),numel(spar));
FAM = FaM(Spar,Sper.^2);
for ip_ = 1:Np
p_ = parray(ip_);
HH = Hp(p_,Spar);
for ij_ = 1:Nj
j_ = jarray(ij_);
LL = Lj(j_,Sper.^2);
HLF = HH.*LL.*FAM;
FF = FF + Napj_(ip_,ij_)*HLF;
end
end
FF = (FF.*conj(FF)); %|f_a|^2
FF = sqrt(FF); %<|f_a|>kx,ky
FF = FF./max(FF(:));
end
\ No newline at end of file
wk/fast_analysis.m
+ 24
− 21
View file @ 968d0648
......@@ -16,9 +16,9 @@ default_plots_options
% resdir ='cheap_CBC_baseline';
% resdir ='test_HEL_closure';
% resdir ='dmax_closure/';
PARTITION = '/misc/gyacomo23_outputs/reduced_fluid_paper/';
% PARTITION = '/misc/gyacomo23_outputs/reduced_fluid_paper/';
% resdir = '/Npol_study/AE_CBC_s0_beta0_P4J2';
resdir = '/Npol_study/RF_CBC_s0_beta0/Npol_11';
% resdir = '/Npol_study/RF_CBC_s0_beta0/Npol_11';
% Triangularity paper
% PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
......@@ -52,10 +52,13 @@ resdir = '/Npol_study/RF_CBC_s0_beta0/Npol_11';
% 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,'/'];
% DATADIR = [PARTITION,resdir,'/'];
% DATADIR = '/home/ahoffman/gyacomo/simulations/ralf/2D_Zpinch_ITG/3x2x64x48x1_no_curvB/';
% DATADIR = '/home/ahoffman/gyacomo/simulations/ralf/3D_Zpinch_ITG/3x2x64x48x16_nocurvB/';
DATADIR = '/home/ahoffman/gyacomo/simulations/ralf/3D_Zpinch_ITG/3x2x64x48x16_nocurvB_-14/';
read_flux_out_XX(DATADIR,1,1);
%%
J0 = 00; J1 = 00;
J0 = 00; J1 = 10;
% Load basic info (grids and time traces)
data = {};
......@@ -91,16 +94,16 @@ if 1
%% 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');
options.TAVG_0 = data.Ts3D(end)/2;
options.TAVG_1 = data.Ts3D(end);
options.TAVG_0 = 250;
options.TAVG_1 = options.TAVG_0+50;
options.NCUT = 5; % Number of cuts for averaging and error estimation
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 = 'n_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'u_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.ST_FIELD = 'T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'n_i T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'Q_{xi}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'G_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
% options.ST_FIELD = 'w_{Ez}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
......@@ -114,7 +117,7 @@ end
if 0
%% 2D field snapshots
% Options
options.INTERP = 1;
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
......@@ -124,23 +127,23 @@ options.TAVG = 1;
% options.NAME = ['N_i^{00}'];
% options.NAME = 'n_e';
% options.NAME = 'upar_i';
% options.NAME = 'n_i';
options.NAME = 'T_i';
% options.NAME = 'Q_{xi}';
% options.NAME = 'v_{Ey}';
% options.NAME = 'w_{Ez}';
% options.NAME = '\omega_z';
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = 'n_i-n_e';
loc =11;
[~,i_] = min(abs(loc - data.grids.y));
options.COMP =i_;
% options.PLAN = '3D';
options.PLAN = '3D';
% options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
options.PLAN = 'kxz'; options.COMP ='avg';
% 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 = [50:500]; options.TAVG = 1;
options.XYZ =[-11 20 -2];
options.TIME = [80 140 190 250]; options.TAVG = 0;
% options.TIME = [50:500]; options.TAVG = 1;
options.RESOLUTION = 256;
fig = photomaton(data,options);
% colormap(gray)
......@@ -281,17 +284,17 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.BWR = 1; % bluewhitered plot or gray
options.CLIMAUTO = 0; % adjust the colormap auto
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = 'w_{Ez}';
% options.NAME = '\psi';
% options.NAME = 'T_i';
options.NAME = 'T_i';
% options.NAME = '\phi^{NZ}';
% options.NAME = ['N_i^{00}'];
% options.NAME = ['N_e^{00}'];
options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
% options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
% options.PLAN = 'xz'; options.COMP ='avg';
% options.PLAN = '3D';
options.XYZ =[-11 20 0];
options.PLAN = '3D';
options.XYZ =[-21 20 0];
options.TIME = data.Ts3D(1:1:end);
% options.TIME = [0:1500];
data.EPS = 0.1;
......
wk/p3_geom_scan_analysis.m
+ 3
− 3
View file @ 968d0648
......@@ -14,15 +14,15 @@ GETFLUXSURFACE = 0;
% partition= '../results/paper_3/';
% Get the scan directory
switch 4
switch 2
case 1 % delta K_T tau=1
casename = 'DIIID rho95 $\tau=1$';
partition= '../results/paper_3/DIIID_tau_1_rho95_geom_scan/';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/';
% % scandir = '3x2x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(2,1)';
% scandir = '3x2x192x48x32_nu_0.1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '3x2x192x48x24_nu_0.1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '3x2x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(2,1)';
scandir = '2_1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '2_1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '5x3x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(4,2)';
scandir = '5x3x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(4,2)';
% scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(2,1)';
......
wk/plot_dist_func_and_trap_pass_el.m
+ 7
− 4
View file @ 968d0648
......@@ -5,10 +5,11 @@ options.SHOW_CURVOP = 0;
[fig, geo_arrays] = plot_metric(data,options);
eps_eff = 1/max((geo_arrays.hatB));
T = [40];
s = linspace(-4,4,64);
x = linspace( 0,4,64);
T = [260];
s = linspace(-3,3,64);
x = linspace( 0,3,64);
[SS,XX,FF,FAM] = compute_fa_2D(data,'i',s,x,T);
% [SS,XX,FF,FAM] = compute_fa_2D_spar_sper(data,'i',s,x,T);
% FF = FF - FAM;
......@@ -20,4 +21,6 @@ contourf(SS,XX,(FF),20)
hold on
% 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
colormap(bluewhitered)
xlabel('$v_\parallel/c_s$');
ylabel('$\mu B_0/T_e$');
\ No newline at end of file
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