From 7be2bab52db05f8d3c04dad1d8930c4b76d8b11e Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 2 Aug 2021 10:34:38 +0200
Subject: [PATCH] Reorganization of the matlab scripts for more readability
---
matlab/COSOlver_matrices_analysis.m | 12 +-
matlab/check_checkpoint.m | 19 -
matlab/compile_results.m | 30 +-
matlab/compute_Sapj.m | 37 --
matlab/conv_thm_2D.m | 14 -
matlab/create_gif.m | 4 +
matlab/create_mov.m | 4 +
matlab/default_plots_options.m | 3 +-
matlab/fort.90 | 6 -
matlab/load_2D_data.m | 6 +-
matlab/load_3D_data.m | 17 +
matlab/load_4D_data.m | 23 +
matlab/load_5D_data.m | 9 +-
matlab/load_grid_data.m | 7 +-
matlab/load_params.m | 12 +-
matlab/load_results.m | 16 +-
matlab/plot_kernels.m | 32 +-
matlab/plots/plot_kperp_spectrum.m | 37 ++
.../plots/plot_radial_transport_and_shear.m | 44 ++
matlab/plots/plot_shear_and_reynold_stress.m | 10 +
matlab/plots/plot_space_time_diagrams.m | 24 ++
matlab/plots/plot_time_evolution_and_gr.m | 64 +++
matlab/post_processing.m | 191 +++++++++
matlab/setup.m | 27 +-
matlab/write_fort90.m | 12 +-
matlab/write_sbash_daint.m | 2 +-
matlab/write_sbash_izar.m | 2 +-
wk/analysis_3D.m | 399 +-----------------
28 files changed, 528 insertions(+), 535 deletions(-)
delete mode 100644 matlab/check_checkpoint.m
delete mode 100644 matlab/compute_Sapj.m
delete mode 100644 matlab/conv_thm_2D.m
delete mode 100644 matlab/fort.90
create mode 100644 matlab/load_3D_data.m
create mode 100644 matlab/load_4D_data.m
create mode 100644 matlab/plots/plot_kperp_spectrum.m
create mode 100644 matlab/plots/plot_radial_transport_and_shear.m
create mode 100644 matlab/plots/plot_shear_and_reynold_stress.m
create mode 100644 matlab/plots/plot_space_time_diagrams.m
create mode 100644 matlab/plots/plot_time_evolution_and_gr.m
create mode 100644 matlab/post_processing.m
diff --git a/matlab/COSOlver_matrices_analysis.m b/matlab/COSOlver_matrices_analysis.m
index 1514d484..083eb5ea 100644
--- a/matlab/COSOlver_matrices_analysis.m
+++ b/matlab/COSOlver_matrices_analysis.m
@@ -1,14 +1,14 @@
addpath(genpath('../matlab')) % ... add
%% Grid configuration
-N = 10; % Frequency gridpoints (Nkr = N/2)
+N = 10; % Frequency gridpoints (Nkx = N/2)
L = 120; % Size of the squared frequency domain
dk = 2*pi/L;
kmax = N/2*dk;
-kr = dk*(0:N/2);
-kz = dk*(0:N/2);
-[KZ, KR]= meshgrid(kz,kr);
-KPERP = sqrt(KR.^2 + KZ.^2);
-kperp = reshape(KPERP,[1,numel(kr)^2]);
+kx = dk*(0:N/2);
+ky = dk*(0:N/2);
+[ky, kx]= meshgrid(ky,kx);
+KPERP = sqrt(kx.^2 + ky.^2);
+kperp = reshape(KPERP,[1,numel(kx)^2]);
kperp = uniquetol(kperp,1e-14);
Nperp = numel(kperp);
COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
diff --git a/matlab/check_checkpoint.m b/matlab/check_checkpoint.m
deleted file mode 100644
index 8265cf4d..00000000
--- a/matlab/check_checkpoint.m
+++ /dev/null
@@ -1,19 +0,0 @@
-function [] = check_checkpoint(cp_ID, OUTPUTS)
-
-cp_fname = OUTPUTS.rstfile0;
-cp_fname = [cp_fname(2:end-1),'_%.2d.h5'];
-cp_fname = sprintf(cp_fname,cp_ID)
-
-tmp = h5read(cp_fname,'/Basic/moments_i');
-Nipj_cp = tmp.real + 1i * tmp.imaginary;
-t_cp = h5readatt(cp_fname,'/Basic','time');
-
-[~,~,Nr,Nz] = size(Nipj_cp);
-
-% ni00_cp = fftshift(ifft2(half_2_full_cc_2D(squeeze(Nipj_cp(1,1,:,:))),'symmetric'));
-ni00_cp = real(fftshift(ifft2(squeeze(Nipj_cp(1,1,:,:)),Nz,Nz)));
-fig = figure;
- pclr=pcolor(transpose(ni00_cp(:,:)));set(pclr, 'edgecolor','none'); colorbar;
- title(['$n_i^{00}$',', $t \approx$', sprintf('%.3d',ceil(t_cp))]);
- xlabel('$i_r$'); ylabel('$i_z$');
-end
\ No newline at end of file
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index eac0efbf..b79637a5 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -1,5 +1,5 @@
CONTINUE = 1;
-JOBNUM = 0; JOBFOUND = 0;
+JOBNUM = JOBNUMMIN; JOBFOUND = 0;
TJOB_SE = []; % Start and end times of jobs
NU_EVOL = []; % evolution of parameter nu between jobs
MU_EVOL = []; % evolution of parameter mu between jobs
@@ -17,7 +17,7 @@ DENS_I_ = [];
TEMP_E_ = [];
TEMP_I_ = [];
Ts0D_ = [];
-Ts2D_ = [];
+Ts3D_ = [];
Ts5D_ = [];
Sipj_ = []; Sepj_ = [];
Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old;
@@ -79,33 +79,33 @@ while(CONTINUE)
end
if W_PHI || W_NA00
- Ts2D_ = cat(1,Ts2D_,Ts2D);
+ Ts3D_ = cat(1,Ts3D_,Ts3D);
end
if W_PHI
- PHI_ = cat(3,PHI_,PHI);
+ PHI_ = cat(4,PHI_,PHI);
end
if W_NA00
- Ni00_ = cat(3,Ni00_,Ni00);
- Ne00_ = cat(3,Ne00_,Ne00);
+ Ni00_ = cat(4,Ni00_,Ni00);
+ Ne00_ = cat(4,Ne00_,Ne00);
end
if W_DENS
- DENS_E_ = cat(3,DENS_E_,DENS_E);
- DENS_I_ = cat(3,DENS_I_,DENS_I);
+ DENS_E_ = cat(4,DENS_E_,DENS_E);
+ DENS_I_ = cat(4,DENS_I_,DENS_I);
end
if W_TEMP
- TEMP_E_ = cat(3,TEMP_E_,TEMP_E);
- TEMP_I_ = cat(3,TEMP_I_,TEMP_I);
+ TEMP_E_ = cat(4,TEMP_E_,TEMP_E);
+ TEMP_I_ = cat(4,TEMP_I_,TEMP_I);
end
if W_NAPJ || W_SAPJ
Ts5D_ = cat(1,Ts5D_,Ts5D);
end
if W_NAPJ
- Nipj_ = cat(5,Nipj_,Nipj);
- Nepj_ = cat(5,Nepj_,Nepj);
+ Nipj_ = cat(6,Nipj_,Nipj);
+ Nepj_ = cat(6,Nepj_,Nepj);
end
if W_SAPJ
- Sipj_ = cat(5,Sipj_,Sipj);
- Sepj_ = cat(5,Sepj_,Sepj);
+ Sipj_ = cat(6,Sipj_,Sipj);
+ Sepj_ = cat(6,Sepj_,Sepj);
end
% Evolution of simulation parameters
@@ -129,7 +129,7 @@ while(CONTINUE)
end
GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_;
Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
-Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_;
+Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts3D = Ts3D_;
DENS_E = DENS_E_; DENS_I = DENS_I_; TEMP_E = TEMP_E_; TEMP_I = TEMP_I_;
clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_
diff --git a/matlab/compute_Sapj.m b/matlab/compute_Sapj.m
deleted file mode 100644
index 0a708aab..00000000
--- a/matlab/compute_Sapj.m
+++ /dev/null
@@ -1,37 +0,0 @@
-function [ Sapj ] = compute_Sapj(p, j, Kr, Kz, Napj, specie, phi, MODEL, GRID)
-%COMPUT_SAPJ compute the non linear term for moment pj specie a
-% perform a convolution by product of inverse fourier transform and then
-% put the results back into k space with an FFT
- Jmax = GRID.jmaxe * (specie=='e')...
- + GRID.jmaxi * (specie=='i');
- %padding
- Pad = 2.0;
- Sapj = zeros(numel(Kr), numel(Kz));
- F = zeros(numel(Kr), numel(Kz));
- G = zeros(numel(Kr), numel(Kz));
- for n = 0:Jmax % Sum over Laguerre
-
- for ikr = 1:numel(Kr)
- for ikz = 1:numel(Kz)
- kr = Kr(ikr); kz = Kz(ikz);
- BA = sqrt(kr^2+kz^2)*...
- (MODEL.sigma_e*sqrt(2) * (specie=='e')...
- + sqrt(2*MODEL.tau_i) * (specie=='i'));
-
- %First conv term
- F(ikr,ikz) = (kz-kr).*phi(ikr,ikz).*kernel(n,BA);
- %Second conv term
- G(ikr,ikz) = 0.0;
- for s = 0:min(n+j,Jmax)
- G(ikr,ikz) = ...
- G(ikr,ikz) + dnjs(n,j,s) .* squeeze(Napj(p+1,s+1,ikr,ikz));
- end
- G(ikr,ikz) = (kz-kr) .* G(ikr,ikz);
- end
- end
- %Conv theorem
- Sapj = Sapj + conv_thm_2D(F,G,Pad);
- end
-end
-
-
diff --git a/matlab/conv_thm_2D.m b/matlab/conv_thm_2D.m
deleted file mode 100644
index 6f2b4452..00000000
--- a/matlab/conv_thm_2D.m
+++ /dev/null
@@ -1,14 +0,0 @@
-function [ conv ] = conv_thm_2D( F, G, Pad )
-%conv_thm_2D computes the convolution between F and G with a zero pading Pad
-% kspace -> real -> product -> kspace
-
-[Nr, Nz] = size(F);
-
-f = ifft2(F,Nr*Pad, Nz*Pad);
-g = ifft2(G,Nr*Pad, Nz*Pad);
-
-conv_pad = fft2(f.*g); % convolution becomes product
-
-conv = conv_pad(1:Nr,1:Nz); % remove padding
-end
-
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index ca3e37b7..cd020f8d 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -2,6 +2,10 @@ title1 = GIFNAME;
FIGDIR = BASIC.RESDIR;
GIFNAME = [FIGDIR, GIFNAME,'.gif'];
+XNAME = latexize(XNAME);
+YNAME = latexize(YNAME);
+FIELDNAME = latexize(FIELDNAME);
+
% Set colormap boundaries
hmax = max(max(max(FIELD)));
hmin = min(min(min(FIELD)));
diff --git a/matlab/create_mov.m b/matlab/create_mov.m
index cebdb625..ca28c3c3 100644
--- a/matlab/create_mov.m
+++ b/matlab/create_mov.m
@@ -1,6 +1,10 @@
title1 = GIFNAME;
FIGDIR = BASIC.RESDIR;
GIFNAME = [FIGDIR, GIFNAME,'.mp4'];
+XNAME = latexize(XNAME);
+YNAME = latexize(YNAME);
+FIELDNAME = latexize(FIELDNAME);
+
vidfile = VideoWriter(GIFNAME,'Uncompressed AVI');
open(vidfile);
% Set colormap boundaries
diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
index 56e500d3..f8d54966 100644
--- a/matlab/default_plots_options.m
+++ b/matlab/default_plots_options.m
@@ -7,4 +7,5 @@ set(0, 'DefaultLineLineWidth', 2.0);
line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];...
[0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]];
line_styles = {'-','-.','--',':'};
-marker_styles = {'^','v','o','s'};
\ No newline at end of file
+marker_styles = {'^','v','o','s'};
+latexize = @(x) ['$',x,'$'];
diff --git a/matlab/fort.90 b/matlab/fort.90
deleted file mode 100644
index 321aad3b..00000000
--- a/matlab/fort.90
+++ /dev/null
@@ -1,6 +0,0 @@
-&GRID
- nr = 64
- Lr = 10
- nz = 64
- Lz = 10
-/
diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m
index 47ef91a6..b013a63d 100644
--- a/matlab/load_2D_data.m
+++ b/matlab/load_2D_data.m
@@ -1,12 +1,12 @@
function [ data, time, dt ] = load_2D_data( filename, variablename )
%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
time = h5read(filename,'/data/var2d/time');
- kr = h5read(filename,'/data/grid/coordkr');
- kz = h5read(filename,'/data/grid/coordkz');
+ kx = h5read(filename,'/data/grid/coordkx');
+ ky = h5read(filename,'/data/grid/coordky');
dt = h5readatt(filename,'/data/input','dt');
cstart= h5readatt(filename,'/data/input','start_iframe2d');
- data = zeros(numel(kr),numel(kz),numel(time));
+ data = zeros(numel(kx),numel(ky),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,it) = tmp.real + 1i * tmp.imaginary;
diff --git a/matlab/load_3D_data.m b/matlab/load_3D_data.m
new file mode 100644
index 00000000..2921908a
--- /dev/null
+++ b/matlab/load_3D_data.m
@@ -0,0 +1,17 @@
+function [ data, time, dt ] = load_3D_data( filename, variablename )
+%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
+ time = h5read(filename,'/data/var3d/time');
+ kx = h5read(filename,'/data/grid/coordkx');
+ ky = h5read(filename,'/data/grid/coordky');
+ z = h5read(filename,'/data/grid/coordz');
+ dt = h5readatt(filename,'/data/input','dt');
+ cstart= h5readatt(filename,'/data/input','start_iframe3d');
+
+ data = zeros(numel(kx),numel(ky),numel(z),numel(time));
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
+
+end
+
diff --git a/matlab/load_4D_data.m b/matlab/load_4D_data.m
new file mode 100644
index 00000000..5b7714e5
--- /dev/null
+++ b/matlab/load_4D_data.m
@@ -0,0 +1,23 @@
+function [ data, time, dt ] = load_4D_data( filename, variablename )
+%LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ
+ time = h5read(filename,'/data/var5d/time');
+ if strcmp(variablename,'moments_e') || strcmp(variablename,'Sepj')
+ p = h5read(filename,'/data/grid/coordp_e');
+ j = h5read(filename,'/data/grid/coordj_e');
+ else
+ p = h5read(filename,'/data/grid/coordp_i');
+ j = h5read(filename,'/data/grid/coordj_i');
+ end
+ kx = h5read(filename,'/data/grid/coordkx');
+ ky = h5read(filename,'/data/grid/coordky');
+
+ dt = h5readatt(filename,'/data/input','dt');
+ cstart= h5readatt(filename,'/data/input','start_iframe5d');
+
+ data = zeros(numel(p),numel(j),numel(kx),numel(ky),numel(time));
+
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
+end
\ No newline at end of file
diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m
index e36d7675..2de05533 100644
--- a/matlab/load_5D_data.m
+++ b/matlab/load_5D_data.m
@@ -8,16 +8,17 @@ function [ data, time, dt ] = load_5D_data( filename, variablename )
p = h5read(filename,'/data/grid/coordp_i');
j = h5read(filename,'/data/grid/coordj_i');
end
- kr = h5read(filename,'/data/grid/coordkr');
- kz = h5read(filename,'/data/grid/coordkz');
+ kx = h5read(filename,'/data/grid/coordkx');
+ ky = h5read(filename,'/data/grid/coordky');
+ z = h5read(filename,'/data/grid/coordz');
dt = h5readatt(filename,'/data/input','dt');
cstart= h5readatt(filename,'/data/input','start_iframe5d');
- data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time));
+ data = zeros(numel(p),numel(j),numel(kx),numel(ky),numel(z),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]);
- data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ data(:,:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end
end
\ No newline at end of file
diff --git a/matlab/load_grid_data.m b/matlab/load_grid_data.m
index 43181577..be964241 100644
--- a/matlab/load_grid_data.m
+++ b/matlab/load_grid_data.m
@@ -1,9 +1,10 @@
-function [ pe, je, pi, ji, kr, kz ] = load_grid_data( filename )
+function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data_3D( filename )
%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ
pe = h5read(filename,'/data/grid/coordp_e');
je = h5read(filename,'/data/grid/coordj_e');
pi = h5read(filename,'/data/grid/coordp_i');
ji = h5read(filename,'/data/grid/coordj_i');
- kr = h5read(filename,'/data/grid/coordkr');
- kz = h5read(filename,'/data/grid/coordkz');
+ kx = h5read(filename,'/data/grid/coordkx');
+ ky = h5read(filename,'/data/grid/coordky');
+ z = h5read(filename,'/data/grid/coordz');
end
\ No newline at end of file
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 28ce1ff9..336c3057 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -8,9 +8,9 @@ PMAXE = h5readatt(filename,'/data/input','pmaxe');
JMAXE = h5readatt(filename,'/data/input','jmaxe');
NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
NU = h5readatt(filename,'/data/input','nu');
-NR = h5readatt(filename,'/data/input','nr');
-NZ = h5readatt(filename,'/data/input','nz');
-L = h5readatt(filename,'/data/input','Lr');
+Nx = h5readatt(filename,'/data/input','Nx');
+Ny = h5readatt(filename,'/data/input','Ny');
+L = h5readatt(filename,'/data/input','Lx');
CLOS = h5readatt(filename,'/data/input','CLOS');
DT_SIM = h5readatt(filename,'/data/input','dt');
MU = h5readatt(filename,'/data/input','mu');
@@ -46,11 +46,11 @@ else
degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
'_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
end
-degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
+degngrad = [degngrad,'_eta_%1.1f_nu_%0.0e_',...
CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[NU,MU]);
+degngrad = sprintf(degngrad,[ETAB/ETAN,NU,MU]);
if ~NON_LIN; degngrad = ['lin_',degngrad]; end
-resolution = [num2str(NR),'x',num2str(NZ/2),'_'];
+resolution = [num2str(Nx),'x',num2str(Ny/2),'_'];
gridname = ['L_',num2str(L),'_'];
PARAMS = [resolution,gridname,degngrad];
% BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
diff --git a/matlab/load_results.m b/matlab/load_results.m
index 412b92a4..fd7e85d6 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -4,7 +4,7 @@ disp(['Loading ',filename])
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input','dt');
-[Pe, Je, Pi, Ji, kr, kz] = load_grid_data(filename);
+[Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi') ,'y');
@@ -21,12 +21,12 @@ if W_GAMMA
end
if W_PHI
- [ PHI, Ts2D, dt2D] = load_2D_data(filename, 'phi');
+ [ PHI, Ts3D, dt3D] = load_3D_data(filename, 'phi');
end
if W_NA00
- [Ni00, Ts2D, dt2D] = load_2D_data(filename, 'Ni00');
- Ne00 = load_2D_data(filename, 'Ne00');
+ [Ni00, Ts3D, dt3D] = load_3D_data(filename, 'Ni00');
+ Ne00 = load_3D_data(filename, 'Ne00');
end
@@ -41,11 +41,11 @@ if W_SAPJ
end
if W_DENS
- [DENS_E, Ts2D, dt2D] = load_2D_data(filename, 'dens_e');
- [DENS_I, Ts2D, dt2D] = load_2D_data(filename, 'dens_i');
+ [DENS_E, Ts3D, dt3D] = load_3D_data(filename, 'dens_e');
+ [DENS_I, Ts3D, dt3D] = load_3D_data(filename, 'dens_i');
end
if W_TEMP
- [TEMP_E, Ts2D, dt2D] = load_2D_data(filename, 'temp_e');
- [TEMP_I, Ts2D, dt2D] = load_2D_data(filename, 'temp_i');
+ [TEMP_E, Ts3D, dt3D] = load_3D_data(filename, 'temp_e');
+ [TEMP_I, Ts3D, dt3D] = load_3D_data(filename, 'temp_i');
end
\ No newline at end of file
diff --git a/matlab/plot_kernels.m b/matlab/plot_kernels.m
index f7729846..f7eede8a 100644
--- a/matlab/plot_kernels.m
+++ b/matlab/plot_kernels.m
@@ -1,16 +1,16 @@
%% Kernels
kmax=7;
nmax=6;
-kr_ = linspace(0,kmax,100);
+kx_ = linspace(0,kmax,100);
figure
for n_ = 0:nmax
- plot(kr_,kernel(n_,kr_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
+ plot(kx_,kernel(n_,kx_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
end
ylim_ = ylim;
-plot(kr_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
-plot(kr_,J0,'-r','DisplayName','$J_0$');
+plot(kx_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+plot(kx_,J0,'-r','DisplayName','$J_0$');
legend('show')
%% Bessels and approx
@@ -22,24 +22,24 @@ figure
for i = 1:4
subplot(2,2,i)
v_ = vperp(i);
- kr_ = linspace(0,kmax,100);
+ kx_ = linspace(0,kmax,100);
- J0 = besselj(0,kr_*v_);
- A1 = 1 - kr_.^2*v_^2/4;
- K1 = zeros(size(kr_));
- K2 = zeros(size(kr_));
+ J0 = besselj(0,kx_*v_);
+ A1 = 1 - kx_.^2*v_^2/4;
+ K1 = zeros(size(kx_));
+ K2 = zeros(size(kx_));
for n_ = 0:nmax1
- K1 = K1 + kernel(n_,kr_).*polyval(LaguerrePoly(n_),v_^2);
+ K1 = K1 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2);
end
for n_ = 0:nmax2
- K2 = K2 + kernel(n_,kr_).*polyval(LaguerrePoly(n_),v_^2);
+ K2 = K2 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2);
end
- plot(kr_,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
- plot(kr_,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
- plot(kr_,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
- plot(kr_,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kx_,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
+ plot(kx_,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
+ plot(kx_,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kx_,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
ylim_ = [-0.5, 1.0];
- plot(kr_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+ plot(kx_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
ylim(ylim_); xlabel('$k$')
legend('show'); grid on; title(['$v = ',num2str(v_),'$'])
end
diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plots/plot_kperp_spectrum.m
new file mode 100644
index 00000000..2370dd02
--- /dev/null
+++ b/matlab/plots/plot_kperp_spectrum.m
@@ -0,0 +1,37 @@
+[~,itstart] = min(abs(Ts3D-tstart));
+[~,itend] = min(abs(Ts3D-tend));
+trange = itstart:itend;
+%full kperp points
+phi_k_2 = reshape(mean(mean(abs(PHI(:,:,:,trange)),3),4).^2,[numel(KX),1]);
+kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
+% interpolated kperps
+nk_noAA = floor(2/3*numel(kx));
+kp_ip = kx;
+[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
+[xn,yn] = pol2cart(thg,rg);
+[ky_s, sortIdx] = sort(ky);
+[xc,yc] = meshgrid(ky_s,kx);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,trange))).^2,3)),xn,yn);
+phi_kp = mean(Z_rth,2);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(Ni00(:,sortIdx,trange))).^2,3)),xn,yn);
+ni00_kp = mean(Z_rth,2);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(Ne00(:,sortIdx,trange))).^2,3)),xn,yn);
+ne00_kp = mean(Z_rth,2);
+%for theorical trends
+a1 = phi_kp(2)*kp_ip(2).^(13/3);
+a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
+fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
+% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
+plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
+plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
+set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on
+xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest')
+title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+' $\mu_{hd}=$',num2str(MU)]});
+xlim([0.1,10]);
+save_figure
+clear Z_rth phi_kp ni_kp Ti_kp
\ No newline at end of file
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
new file mode 100644
index 00000000..ba865250
--- /dev/null
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -0,0 +1,44 @@
+%Compute steady radial transport
+tend = Ts0D(end); tstart = tend - TAVG;
+[~,its0D] = min(abs(Ts0D-tstart));
+[~,ite0D] = min(abs(Ts0D-tend));
+SCALE = (2*pi/Nx/Ny)^2;
+gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
+gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
+% Compute steady shearing rate
+tend = Ts3D(end); tstart = tend - TAVG;
+[~,its2D] = min(abs(Ts3D-tstart));
+[~,ite2D] = min(abs(Ts3D-tend));
+shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1));
+shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1));
+% plots
+fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
+ subplot(311)
+% yyaxis left
+ plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
+ plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
+ 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
+ grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$')
+ ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+ %
+ subplot(312)
+ clr = line_colors(1,:);
+ lstyle = line_styles(1);
+% plt = @(x_) mean(x_,1);
+ plt = @(x_) x_(1,:);
+ plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
+ plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
+ plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
+ plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on;
+ plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
+ 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
+ ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
+ grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
+ subplot(313)
+ [TY,TX] = meshgrid(x,Ts3D);
+ pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
+ caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
+save_figure
\ No newline at end of file
diff --git a/matlab/plots/plot_shear_and_reynold_stress.m b/matlab/plots/plot_shear_and_reynold_stress.m
new file mode 100644
index 00000000..e951208a
--- /dev/null
+++ b/matlab/plots/plot_shear_and_reynold_stress.m
@@ -0,0 +1,10 @@
+% trange = 100:200;
+figure;
+plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3))));
+plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on;
+plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on;
+% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on;
+plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on;
+plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on;
+% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on;
+xlim([-L/2,L/2]); xlabel('$x/\rho_s$'); grid on; legend('show')
\ No newline at end of file
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m
new file mode 100644
index 00000000..3e263a0a
--- /dev/null
+++ b/matlab/plots/plot_space_time_diagrams.m
@@ -0,0 +1,24 @@
+[~,itstart] = min(abs(Ts3D-tstart));
+[~,itend] = min(abs(Ts3D-tend));
+trange = itstart:itend;
+[TY,TX] = meshgrid(x,Ts3D(trange));
+fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
+ subplot(211)
+% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+ pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+ shading interp
+ colormap hot;
+ caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
+ caxis([0.0,cmax]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
+ xticks([]); ylabel('$x/\rho_s$')
+% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
+ title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+ subplot(212)
+ pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+ fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2)));
+ caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
+ xlabel('$t c_s/R$'), ylabel('$x/\rho_s$')
+% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
+save_figure
\ No newline at end of file
diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plots/plot_time_evolution_and_gr.m
new file mode 100644
index 00000000..e016520b
--- /dev/null
+++ b/matlab/plots/plot_time_evolution_and_gr.m
@@ -0,0 +1,64 @@
+fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
+set(gcf, 'Position', [100, 100, 900, 800])
+subplot(111);
+ suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+ subplot(421);
+ for ip = 1:Pe_max
+ for ij = 1:Je_max
+ plt = @(x) squeeze(x(ip,ij,:));
+ plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
+ clr = line_colors(min(ip,numel(line_colors(:,1))),:);
+ lstyle = line_styles(min(ij,numel(line_styles)));
+ semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,...
+ 'Color',clr,'LineStyle',lstyle{1}); hold on;
+ end
+ end
+ grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
+ subplot(423)
+ for ip = 1:Pi_max
+ for ij = 1:Ji_max
+ plt = @(x) squeeze(x(ip,ij,:));
+ plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
+ clr = line_colors(min(ip,numel(line_colors(:,1))),:);
+ lstyle = line_styles(min(ij,numel(line_styles)));
+ semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,...
+ 'Color',clr,'LineStyle',lstyle{1}); hold on;
+ end
+ end
+ grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
+ subplot(222)
+ plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on;
+ plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2);
+ legend(['Gyro. flux';'Part. flux']);
+ grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
+ if(~isnan(max(max(g_I(1,:,:)))))
+ subplot(223)
+ plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on;
+ plot(kx,max(g_II(:,1,:),[],3),'x-','DisplayName','Second. instability'); hold on;
+ plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
+ grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
+ ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]);
+ shearplot = 426; phiplot = 428;
+ else
+ shearplot = 223; phiplot = 224;
+ end
+ subplot(shearplot)
+ plt = @(x) mean(x,1);
+ clr = line_colors(min(ip,numel(line_colors(:,1))),:);
+ lstyle = line_styles(min(ij,numel(line_styles)));
+ plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s)$'); hold on;
+ plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on;
+ plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on;
+ plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on;
+ grid on; xlabel('$t c_s/R$'); ylabel('$shear$');
+ subplot(phiplot)
+ clr = line_colors(min(ip,numel(line_colors(:,1))),:);
+ lstyle = line_styles(min(ij,numel(line_styles)));
+ plot(Ts3D,plt(phi_maxx_maxy),'DisplayName','$\max_{r,z}(\phi)$'); hold on;
+ plot(Ts3D,plt(phi_maxx_avg),'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on;
+ plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
+ plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
+ grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$');
+save_figure
\ No newline at end of file
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
new file mode 100644
index 00000000..ea25ef89
--- /dev/null
+++ b/matlab/post_processing.m
@@ -0,0 +1,191 @@
+
+%% Retrieving max polynomial degree and sampling info
+Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
+Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
+Ns5D = numel(Ts5D);
+Ns3D = numel(Ts3D);
+% renaming and reshaping quantity of interest
+Ts5D = Ts5D';
+Ts3D = Ts3D';
+
+%% Build grids
+Nkx = numel(kx); Nky = numel(ky);
+[KY,KX] = meshgrid(ky,kx);
+Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
+dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
+KPERP2 = KY.^2+KX.^2;
+[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
+
+Lk = max(Lkx,Lky);
+Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z);
+dx = 2*pi/Lk; dy = 2*pi/Lk; dz = q0*2*pi/Nz;
+x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
+y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
+z = dz * (1:Nz);
+[Y_XY,X_XY] = meshgrid(y,x);
+[Z_XZ,X_XZ] = meshgrid(z,x);
+[Z_YZ,Y_YZ] = meshgrid(z,y);
+
+%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+disp('Analysis :')
+disp('- iFFT')
+% IFFT (Lower case = real space, upper case = frequency space)
+ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density
+ni00 = zeros(Nx,Ny,Nz,Ns3D);
+dzTe = zeros(Nx,Ny,Nz,Ns3D);
+dzTi = zeros(Nx,Ny,Nz,Ns3D);
+dzni = zeros(Nx,Ny,Nz,Ns3D);
+np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density
+si00 = zeros(Nx,Ny,Nz,Ns5D);
+phi = zeros(Nx,Ny,Nz,Ns3D);
+dens_e = zeros(Nx,Ny,Nz,Ns3D);
+dens_i = zeros(Nx,Ny,Nz,Ns3D);
+temp_e = zeros(Nx,Ny,Nz,Ns3D);
+temp_i = zeros(Nx,Ny,Nz,Ns3D);
+drphi = zeros(Nx,Ny,Nz,Ns3D);
+dzphi = zeros(Nx,Ny,Nz,Ns3D);
+dr2phi = zeros(Nx,Ny,Nz,Ns3D);
+
+for it = 1:numel(Ts3D)
+ for iz = 1:numel(z)
+ NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
+ ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
+ ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
+ phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
+ drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny)));
+ dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
+ dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
+ if(W_DENS && W_TEMP)
+ DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it);
+ TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it);
+ dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
+ dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
+ dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
+ dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
+ dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
+ temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
+ temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
+ end
+ end
+end
+
+Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space
+
+for it = 1:numel(Ts5D)
+ [~, it2D] = min(abs(Ts3D-Ts5D(it)));
+ Np_i(:,:,it) = 0;
+ for ij = 1:Nji
+ Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
+ Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it));
+ end
+
+ np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny)));
+end
+
+% Post processing
+disp('- post processing')
+% gyrocenter and particle flux from real space
+GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni drphi>
+GFlux_yi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni dzphi>
+GFlux_xe = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne drphi>
+GFlux_ye = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne dzphi>
+% Hermite energy spectrum
+epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D);
+epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D);
+
+phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
+phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
+phi_maxx_avg = zeros(Nz,Ns3D); % Time evol. of the norm of phi
+phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
+
+shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
+shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
+shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
+shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
+
+Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj
+Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % .
+
+% Kperp spectrum interpolation
+%full kperp points
+kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
+% interpolated kperps
+nk_noAA = floor(2/3*numel(kx));
+kp_ip = kx;
+[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
+[xn,yn] = pol2cart(thg,rg);
+[ky_s, sortIdx] = sort(ky);
+[xc,yc] = meshgrid(ky_s,kx);
+phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D);
+%
+for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
+ for iz = 1:numel(z)
+ NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
+ phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it))));
+ phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it))));
+ phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it))));
+ phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
+
+ shear_maxx_maxy(iz,it) = max( max(squeeze(-(dr2phi(:,:,iz,it)))));
+ shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dr2phi(:,:,iz,it)))));
+ shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dr2phi(:,:,iz,it)))));
+ shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dr2phi(:,:,iz,it)))));
+
+ GFlux_xi(iz,it) = sum(sum(ni00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
+ GFlux_yi(iz,it) = sum(sum(-ni00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
+ GFlux_xe(iz,it) = sum(sum(ne00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
+ GFlux_ye(iz,it) = sum(sum(-ne00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
+
+ Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn);
+ phi_kp_t(:,iz,it) = mean(Z_rth,2);
+ end
+end
+%
+for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
+ [~, it2D] = min(abs(Ts3D-Ts5D(it)));
+ Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
+ Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
+ epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4);
+ epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4);
+end
+
+%% Compute primary instability growth rate
+disp('- growth rate')
+% Find max value of transport (end of linear mode)
+[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
+[~,itmax] = min(abs(Ts3D-tmax));
+tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax);
+[~,its3D_lin] = min(abs(Ts3D-tstart));
+[~,ite3D_lin] = min(abs(Ts3D-tend));
+
+g_I = zeros(Nkx,Nky,Nz);
+for ikx = 1:Nkx
+ for iky = 1:Nky
+ for iz = 1:Nz
+ [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+ end
+ end
+end
+[gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
+kmax_I = abs(ky(ikmax_I));
+Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2;
+
+%% Compute secondary instability growth rate
+disp('- growth rate')
+% Find max value of transport (end of linear mode)
+% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
+% [~,itmax] = min(abs(Ts2D-tmax));
+% tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
+[~,its3D_lin] = min(abs(Ts3D-tstart));
+[~,ite3D_lin] = min(abs(Ts3D-tend));
+
+g_II = zeros(Nkx,Nky);
+for ikx = 1:Nkx
+ for iky = 1
+ for iz = 1:Nz
+ [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+ end
+ end
+end
+[gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
+kmax_II = abs(kx(ikmax_II));
diff --git a/matlab/setup.m b/matlab/setup.m
index 92352d9d..517c9217 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -5,11 +5,12 @@ GRID.pmaxe = PMAXE; % Electron Hermite moments
GRID.jmaxe = JMAXE; % Electron Laguerre moments
GRID.pmaxi = PMAXI; % Ion Hermite moments
GRID.jmaxi = JMAXI; % Ion Laguerre moments
-GRID.Nr = N; % r grid resolution
-GRID.Lr = L; % r length
-GRID.Nz = N * (1-KREQ0) + KREQ0; % z ''
-GRID.Lz = L * (1-KREQ0); % z ''
-GRID.kpar = KPAR;
+GRID.Nx = N; % x grid resolution
+GRID.Lx = L; % x length
+GRID.Ny = N * (1-KXEQ0) + KXEQ0; % y ''
+GRID.Ly = L * (1-KXEQ0); % y ''
+GRID.Nz = Nz; % z resolution
+GRID.q0 = q0; % q factor
% Model parameters
MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
@@ -35,11 +36,13 @@ MODEL.eta_n = ETAN; % source term kappa for HW
MODEL.eta_T = ETAT; % Temperature
MODEL.eta_B = ETAB; % Magnetic
MODEL.lambdaD = LAMBDAD;
-% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
+% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KX0KH),'_A_',num2str(A0KH)]; end;
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
INITIAL.INIT_ZF = INIT_ZF;
+if (WIPE_TURB); INITIAL.wipe_turb = '.true.'; else; INITIAL.wipe_turb = '.false.';end;
+if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end;
INITIAL.init_background = (INIT_ZF>0)*ZF_AMP;
INITIAL.init_noiselvl = NOISE0;
INITIAL.iseed = 42;
@@ -74,12 +77,17 @@ else
'_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
end
degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
- CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
+ CONAME,'_mu_%0.0e'];
degngrad = sprintf(degngrad,[NU,MU]);
if ~NON_LIN; degngrad = ['lin_',degngrad]; end
-resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
-gridname = ['L_',num2str(L),'_'];
+if (Nz == 1)
+ resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'_'];
+ gridname = ['L_',num2str(L),'_'];
+else
+ resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'x',num2str(GRID.Nz),'_'];
+ gridname = ['L_',num2str(L),'_q_',num2str(q),'_'];
+end
if (exist('PREFIX','var') == 0); PREFIX = []; end;
if (exist('SUFFIX','var') == 0); SUFFIX = []; end;
PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX];
@@ -98,6 +106,7 @@ if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT);
OUTPUTS.nsave_1d = -1;
OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
+OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT);
OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);
OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT);
if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 19ccdf5d..366ee163 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -16,17 +16,19 @@ fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),'\n']);
fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']);
fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']);
fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']);
-fprintf(fid,[' Nr = ', num2str(GRID.Nr),'\n']);
-fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']);
+fprintf(fid,[' Nx = ', num2str(GRID.Nx),'\n']);
+fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']);
+fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']);
+fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
-fprintf(fid,[' Lz = ', num2str(GRID.Lz),'\n']);
-fprintf(fid,[' kpar = ', num2str(GRID.kpar),'\n']);
+fprintf(fid,[' q0 = ', num2str(GRID.q0),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
+fprintf(fid,[' nsave_3d = ', num2str(OUTPUTS.nsave_3d),'\n']);
fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
fprintf(fid,[' nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']);
fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
@@ -67,6 +69,8 @@ fprintf(fid,'/\n');
fprintf(fid,'&INITIAL_CON\n');
fprintf(fid,[' INIT_NOISY_PHI =', INITIAL.init_noisy_phi,'\n']);
fprintf(fid,[' INIT_ZF =', num2str(INITIAL.INIT_ZF),'\n']);
+fprintf(fid,[' WIPE_TURB =', INITIAL.wipe_turb,'\n']);
+fprintf(fid,[' INIT_BLOB =', INITIAL.init_blob,'\n']);
fprintf(fid,[' init_background =', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed =', num2str(INITIAL.iseed),'\n']);
diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m
index 98c1e359..62a48ac2 100644
--- a/matlab/write_sbash_daint.m
+++ b/matlab/write_sbash_daint.m
@@ -48,7 +48,7 @@ fprintf(fid,[...
'module load cray-mpich\n',...
'module load craype-x86-skylake\n',...
'module load cray-fftw\n',...
-'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]);
+'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]);
%'srun ./../../../bin/helaz']);
fclose(fid);
diff --git a/matlab/write_sbash_izar.m b/matlab/write_sbash_izar.m
index 01c7334c..02224f47 100644
--- a/matlab/write_sbash_izar.m
+++ b/matlab/write_sbash_izar.m
@@ -38,7 +38,7 @@ fprintf(fid,[...
'module load intel-mpi/2019.5.281\n',...
'module load fftw/3.3.8-mpi-openmp\n',...
'module load hdf5/1.10.6-mpi\n',...
-'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]);
+'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]);
fclose(fid);
system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 6fc1d667..397c03bc 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,12 +1,8 @@
addpath(genpath('../matlab')) % ... add
-for i_ = 1
-% for ETA_ =[0.6:0.1:0.9]
-%% Load results
+addpath(genpath('../matlab/plots')) % ... add
+%% Directory of the simulation
if 1% Local results
-outfile ='';
-outfile ='Blob_diffusion/100x50_L_60_P_2_J_1_eta_Inf_nu_1e-01_DGGK_mu_0e+00';
-% outfile ='test_3D/100x50x10_L_60_q0_1_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_2e-03';
-% outfile ='test_3D/100x50_L_60_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_0e+00';
+outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
@@ -23,200 +19,13 @@ BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
-%%
+%% Load the results
+% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 20;
-compile_results_3D %Compile the results from first output found to JOBNUMMAX if existing
-
-%% Retrieving max polynomial degree and sampling info
-Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
-Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
-Ns5D = numel(Ts5D);
-Ns3D = numel(Ts3D);
-% renaming and reshaping quantity of interest
-Ts5D = Ts5D';
-Ts3D = Ts3D';
-
-%% Build grids
-Nkx = numel(kx); Nky = numel(ky);
-[KY,KX] = meshgrid(ky,kx);
-Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
-dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
-KPERP2 = KY.^2+KX.^2;
-[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
-
-Lk = max(Lkx,Lky);
-Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z);
-dx = 2*pi/Lk; dy = 2*pi/Lk; dz = q0*2*pi/Nz;
-x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
-y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
-z = dz * (1:Nz);
-[Y_XY,X_XY] = meshgrid(y,x);
-[Z_XZ,X_XZ] = meshgrid(z,x);
-[Z_YZ,Y_YZ] = meshgrid(z,y);
-
-%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-disp('Analysis :')
-disp('- iFFT')
-% IFFT (Lower case = real space, upper case = frequency space)
-ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density
-ni00 = zeros(Nx,Ny,Nz,Ns3D);
-dzTe = zeros(Nx,Ny,Nz,Ns3D);
-dzTi = zeros(Nx,Ny,Nz,Ns3D);
-dzni = zeros(Nx,Ny,Nz,Ns3D);
-np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density
-si00 = zeros(Nx,Ny,Nz,Ns5D);
-phi = zeros(Nx,Ny,Nz,Ns3D);
-dens_e = zeros(Nx,Ny,Nz,Ns3D);
-dens_i = zeros(Nx,Ny,Nz,Ns3D);
-temp_e = zeros(Nx,Ny,Nz,Ns3D);
-temp_i = zeros(Nx,Ny,Nz,Ns3D);
-drphi = zeros(Nx,Ny,Nz,Ns3D);
-dzphi = zeros(Nx,Ny,Nz,Ns3D);
-dr2phi = zeros(Nx,Ny,Nz,Ns3D);
-
-for it = 1:numel(Ts3D)
- for iz = 1:numel(z)
- NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
- ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
- ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
- phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
- drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny)));
- dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
- dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
- if(W_DENS && W_TEMP)
- DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it);
- TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it);
- dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
- dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
- dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
- dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
- dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
- temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
- temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
- end
- end
-end
-
-Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space
+compile_results %Compile the results from first output found to JOBNUMMAX if existing
-for it = 1:numel(Ts5D)
- [~, it2D] = min(abs(Ts3D-Ts5D(it)));
- Np_i(:,:,it) = 0;
- for ij = 1:Nji
- Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
- Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it));
- end
-
- np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny)));
-end
-
-% Post processing
-disp('- post processing')
-% gyrocenter and particle flux from real space
-GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni drphi>
-GFlux_yi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni dzphi>
-GFlux_xe = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne drphi>
-GFlux_ye = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne dzphi>
-% Hermite energy spectrum
-epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D);
-epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D);
-
-phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_maxx_avg = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-
-shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-
-Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj
-Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % .
-
-% Kperp spectrum interpolation
-%full kperp points
-kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
-% interpolated kperps
-nk_noAA = floor(2/3*numel(kx));
-kp_ip = kx;
-[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
-[xn,yn] = pol2cart(thg,rg);
-[ky_s, sortIdx] = sort(ky);
-[xc,yc] = meshgrid(ky_s,kx);
-phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D);
-%
-for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
- for iz = 1:numel(z)
- NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
- phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it))));
- phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it))));
- phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it))));
- phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
-
- shear_maxx_maxy(iz,it) = max( max(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dr2phi(:,:,iz,it)))));
-
- GFlux_xi(iz,it) = sum(sum(ni00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_yi(iz,it) = sum(sum(-ni00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_xe(iz,it) = sum(sum(ne00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_ye(iz,it) = sum(sum(-ne00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
-
- Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn);
- phi_kp_t(:,iz,it) = mean(Z_rth,2);
- end
-end
-%
-for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
- [~, it2D] = min(abs(Ts3D-Ts5D(it)));
- Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
- Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
- epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4);
- epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4);
-end
-
-%% Compute primary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-[~,itmax] = min(abs(Ts3D-tmax));
-tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax);
-[~,its3D_lin] = min(abs(Ts3D-tstart));
-[~,ite3D_lin] = min(abs(Ts3D-tend));
-
-g_I = zeros(Nkx,Nky,Nz);
-for ikx = 1:Nkx
- for iky = 1:Nky
- for iz = 1:Nz
- [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
- end
- end
-end
-[gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
-kmax_I = abs(ky(ikmax_I));
-Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2;
-
-%% Compute secondary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-% [~,itmax] = min(abs(Ts2D-tmax));
-% tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
-[~,its3D_lin] = min(abs(Ts3D-tstart));
-[~,ite3D_lin] = min(abs(Ts3D-tend));
-
-g_II = zeros(Nkx,Nky);
-for ikx = 1:Nkx
- for iky = 1
- for iz = 1:Nz
- [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
- end
- end
-end
-[gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
-kmax_II = abs(kx(ikmax_II));
+%% Post-processing
+post_processing
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
@@ -225,205 +34,32 @@ FMT = '.fig';
if 1
%% Time evolutions and growth rate
-fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
-set(gcf, 'Position', [100, 100, 900, 800])
-subplot(111);
- suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
- subplot(421);
- for ip = 1:Pe_max
- for ij = 1:Je_max
- plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,...
- 'Color',clr,'LineStyle',lstyle{1}); hold on;
- end
- end
- grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
- subplot(423)
- for ip = 1:Pi_max
- for ij = 1:Ji_max
- plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,...
- 'Color',clr,'LineStyle',lstyle{1}); hold on;
- end
- end
- grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
- subplot(222)
- plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on;
- plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2);
- legend(['Gyro. flux';'Part. flux']);
- grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
- if(~isnan(max(max(g_I(1,:,:)))))
- subplot(223)
- plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on;
- plot(kx,max(g_II(:,1,:),[],3),'x-','DisplayName','Second. instability'); hold on;
- plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
- grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
- ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]);
- shearplot = 426; phiplot = 428;
- else
- shearplot = 223; phiplot = 224;
- end
- subplot(shearplot)
- plt = @(x) mean(x,1);
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s)$'); hold on;
- plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on;
- plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on;
- plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$shear$');
- subplot(phiplot)
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- plot(Ts3D,plt(phi_maxx_maxy),'DisplayName','$\max_{r,z}(\phi)$'); hold on;
- plot(Ts3D,plt(phi_maxx_avg),'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on;
- plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
- plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$');
-save_figure
+plot_time_evolution_and_gr
end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
TAVG = 5000; % Averaging time duration
-%Compute steady radial transport
-tend = Ts0D(end); tstart = tend - TAVG;
-[~,its0D] = min(abs(Ts0D-tstart));
-[~,ite0D] = min(abs(Ts0D-tend));
-SCALE = (2*pi/Nx/Ny)^2;
-gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
-gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
-% Compute steady shearing rate
-tend = Ts3D(end); tstart = tend - TAVG;
-[~,its2D] = min(abs(Ts3D-tstart));
-[~,ite2D] = min(abs(Ts3D-tend));
-shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1));
-shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1));
-% plots
-fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
- subplot(311)
-% yyaxis left
- plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
- plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
- 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
- grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$')
- ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
- %
- subplot(312)
- clr = line_colors(1,:);
- lstyle = line_styles(1);
-% plt = @(x_) mean(x_,1);
- plt = @(x_) x_(1,:);
- plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
- plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
- plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
- plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on;
- plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
- 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
- ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
- grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
- subplot(313)
- [TY,TX] = meshgrid(x,Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
- caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
-save_figure
+plot_radial_transport_and_shear
end
if 0
%% Space time diagramms
-tstart = 0; tend = Ts3D(end);
-[~,itstart] = min(abs(Ts3D-tstart));
-[~,itend] = min(abs(Ts3D-tend));
-trange = itstart:itend;
-[TY,TX] = meshgrid(x,Ts3D(trange));
-fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
- subplot(211)
-% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- shading interp
- colormap hot;
- caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
- caxis([0.0,0.01]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
- xticks([]); ylabel('$x/\rho_s$')
-% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
- subplot(212)
- pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2)));
- caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
- xlabel('$t c_s/R$'), ylabel('$x/\rho_s$')
-% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
-save_figure
+cmax = 0.01 % max of the colorbar for transport
+tstart = 0; tend = Ts3D(end); % time window
+plot_space_time_diagrams
end
if 0
%% Averaged shear and Reynold stress profiles
trange = its2D:ite2D;
-% trange = 100:200;
-figure;
-plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3))));
-plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on;
-plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on;
-% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on;
-plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on;
-plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on;
-% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on;
-xlim([-L/2,L/2]); xlabel('$x/\rho_s$'); grid on; legend('show')
+plot_shear_and_reynold_stress
end
if 0
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 0.8*Ts3D(end); tend = Ts3D(end);
-[~,itstart] = min(abs(Ts3D-tstart));
-[~,itend] = min(abs(Ts3D-tend));
-trange = itstart:itend;
-%full kperp points
-phi_k_2 = reshape(mean(mean(abs(PHI(:,:,:,trange)),3),4).^2,[numel(KX),1]);
-kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
-% interpolated kperps
-nk_noAA = floor(2/3*numel(kx));
-kp_ip = kx;
-[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
-[xn,yn] = pol2cart(thg,rg);
-[ky_s, sortIdx] = sort(ky);
-[xc,yc] = meshgrid(ky_s,kx);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,trange))).^2,3)),xn,yn);
-phi_kp = mean(Z_rth,2);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(Ni00(:,sortIdx,trange))).^2,3)),xn,yn);
-ni00_kp = mean(Z_rth,2);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(Ne00(:,sortIdx,trange))).^2,3)),xn,yn);
-ne00_kp = mean(Z_rth,2);
-%for theorical trends
-a1 = phi_kp(2)*kp_ip(2).^(13/3);
-a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
-fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
-% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
-plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
-plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
-set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on
-xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest')
-title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
-', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-' $\mu_{hd}=$',num2str(MU)]});
-xlim([0.1,10]);
-save_figure
-clear Z_rth phi_kp ni_kp Ti_kp
+tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
+plot_kperp_spectrum
end
if 0
@@ -506,5 +142,4 @@ save_figure
end
%%
-% ZF_fourier_analysis
-end
\ No newline at end of file
+% ZF_fourier_analysis
\ No newline at end of file
--
GitLab