diff --git a/matlab/compute/compute_polarisation.m b/matlab/compute/compute_polarisation.m
new file mode 100644
index 0000000000000000000000000000000000000000..c1ac9f5bb18c900edf335fb029cbca1409fbfdfb
--- /dev/null
+++ b/matlab/compute/compute_polarisation.m
@@ -0,0 +1,32 @@
+function [ pola ] = compute_polarisation( data )
+%compute_polarisation computes the polarisation term (1-Gamma0)phi with
+%ordering up to Jmax
+% Detailed explanation goes here
+PHI = data.PHI;
+T = data.Ts3D;
+Ji = data.Jmaxi;
+pola = zeros(size(PHI));
+
+KN_ = zeros(data.Nky,data.Nkx,data.Nz);
+
+[KX,KY] = meshgrid(data.kx,data.ky);
+KP = sqrt(KX.^2+KY.^2);
+
+GAMMA2_ = 0.*KN_;
+
+for in = 1:Ji+1
+ for iz = 1:data.Nz
+ GAMMA2_(:,:,iz) = GAMMA2_(:,:,iz) + kernel(in-1,KP*sqrt(2)).^2;
+ end
+end
+
+
+for it = 1:numel(T)
+
+ pola(:,:,:,it) = (1-GAMMA2_).*PHI(:,:,:,it);
+
+end
+
+
+end
+
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 2e03ee8e80e1c99a35b7e571cf5c74564151605c..61256c51ccbf9a3c5f78571d60f487d3ddfe3c8a 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,9 +1,9 @@
% Options
SHOW_FILM = 0;
-U_TIME =2400; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME = 50; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin = 1000;
-dt_ = 0.2;
+Tfin = 50;
+dt_ = 0.05;
Nstep = ceil(Tfin/dt_);
% Init tracers
Np = 100; %number of tracers
@@ -48,7 +48,7 @@ end
% position grid and velocity field
[YY_, XX_ ,ZZ_] = meshgrid(data.y,data.x,data.z);
-[KY,KX] = meshgrid(data.ky,data.kx);
+[KX,KY] = meshgrid(data.kx,data.ky);
Ux = zeros(size(XX_));
Uy = zeros(size(XX_));
Uz = zeros(size(XX_));
@@ -57,9 +57,12 @@ ni = zeros(size(XX_));
[~,itu_] = min(abs(U_TIME-data.Ts3D));
% computing the frozen velocity field
for iz = 1:data.Nz
- Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
- Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
- ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
+% Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+% Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+% ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
+ Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
+ Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
+ ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))'*sqrt(data.scale);
end
%% FILM options
@@ -89,9 +92,9 @@ while(t_<Tfin && it <= Nstep)
if Evolve_U && (itu_old ~= itu_)
% updating the velocity field and density field
for iz = 1:data.Nz
- Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
- Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
- ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
+ Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
+ Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))'*sqrt(data.scale);
+ ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))'*sqrt(data.scale);
end
end
% evolve each tracer
@@ -207,7 +210,7 @@ while(t_<Tfin && it <= Nstep)
if SHOW_FILM && (~Evolve_U || (itu_old ~= itu_))
% updating the velocity field
clf(fig);
- F2P = real(ifft2(data.PHI(:,:,iz,itu_),data.Nx,data.Ny));
+ F2P = ifourier_GENE(data.PHI(:,:,iz,itu_))';
scale = max(max(abs(F2P))); % Scaling to normalize
pclr = pcolor(XX_,YY_,F2P/scale);
colormap(bluewhitered);
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 2fdb98a43f7ba94286793fe048cd9a8701b45225..f3ae8237a4433fd3e087fd72e6a27d16fa65ed71 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -3,22 +3,22 @@ function [ DATA ] = load_gene_data( folder )
namelist = read_namelist([folder,'parameters.dat']);
DATA.namelist = namelist;
%% Grid
-file = 'coord.dat.h5';
-DATA.vp = h5read([folder,file],'/coord/vp');
+coofile = 'coord.dat.h5';
+DATA.vp = h5read([folder,coofile],'/coord/vp');
DATA.Nvp = numel(DATA.vp);
-DATA.mu = h5read([folder,file],'/coord/mu');
+DATA.mu = h5read([folder,coofile],'/coord/mu');
DATA.Nmu = numel(DATA.mu);
-DATA.kx = h5read([folder,file],'/coord/kx');
+DATA.kx = h5read([folder,coofile],'/coord/kx');
DATA.Nkx = numel(DATA.kx);
DATA.Nx = DATA.Nkx;
-DATA.ky = h5read([folder,file],'/coord/ky');
+DATA.ky = h5read([folder,coofile],'/coord/ky');
DATA.Nky = numel(DATA.ky);
DATA.Ny = DATA.Nky*2-1;
-DATA.z = h5read([folder,file],'/coord/z');
+DATA.z = h5read([folder,coofile],'/coord/z');
DATA.Nz = numel(DATA.z);
dkx = DATA.kx(2); dky = DATA.ky(2);
@@ -27,33 +27,37 @@ x = linspace(-Lx/2,Lx/2,DATA.Nx+1); x = x(1:end-1);
y = linspace(-Ly/2,Ly/2,DATA.Ny+1); y = y(1:end-1);
DATA.x = x; DATA.y = y; DATA.Lx = Lx; DATA.Ly = Ly;
%% Transport
-file = 'nrg.dat.h5';
-DATA.Ts0D = h5read([folder,file],'/nrgions/time');
-DATA.PGAMMA_RI = h5read([folder,file],'/nrgions/Gamma_es');
-DATA.HFLUX_X = h5read([folder,file],'/nrgions/Q_es');
+nrgfile = 'nrg.dat.h5';
+% nrgfile = 'nrg_1.h5';
+DATA.Ts0D = h5read([folder,nrgfile],'/nrgions/time');
+DATA.PGAMMA_RI = h5read([folder,nrgfile],'/nrgions/Gamma_es');
+DATA.HFLUX_X = h5read([folder,nrgfile],'/nrgions/Q_es');
%% fields and moments
-file = 'mom_ions.dat.h5';
-DATA.Ts3D = h5read([folder,file],'/mom_ions/time');
+phifile = 'field.dat.h5';
+% phifile = 'field_1.h5';
+DATA.Ts3D = h5read([folder,phifile],'/field/time');
DATA.DENS_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
DATA.TPER_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
DATA.TPAR_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
DATA.PHI = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+momfile = 'mom_ions.dat.h5';
+% momfile = 'mom_ions_1.h5';
for jt = 1:numel(DATA.Ts3D)
t = DATA.Ts3D(jt);
[~, it] = min(abs(DATA.Ts3D-t));
- tmp = h5read([folder,file],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+ tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
- tmp = h5read([folder,file],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
+ tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
- tmp = h5read([folder,file],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
+ tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
- tmp = h5read([folder,'field.dat.h5'],['/field/phi/',sprintf('%10.10d',it-1)]);
+ tmp = h5read([folder,phifile],['/field/phi/',sprintf('%10.10d',it-1)]);
DATA.PHI(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
end
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index a56fc536a9082d3f411f463115f57cd15e5af8c5..7ad93c0b6bbb408d6f7420c8312e990407d0b3af 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -68,7 +68,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
ylabel('$Q_x$')
-% ylim([0,2*abs(Qx_infty_avg)]);
+ ylim([0,5*abs(Qx_infty_avg)]);
xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
diff --git a/matlab/setup.m b/matlab/setup.m
index 41d31c2112651cc91b02239a586e584f91504fd2..6288c1de688fdc68d47fe40347b7544c096bf008 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -25,6 +25,7 @@ MODEL.KIN_E = KIN_E;
if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
+MODEL.N_HD = N_HD;
MODEL.mu_z = MU_Z;
MODEL.mu_p = MU_P;
MODEL.mu_j = MU_J;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 26439459c73f738b6faf7c9f9e2e363be52108f1..17838bc320df017124f422c0407357e41b41e1ee 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -57,6 +57,7 @@ fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']);
fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']);
fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']);
fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']);
+fprintf(fid,[' N_HD = ', num2str(MODEL.N_HD),'\n']);
fprintf(fid,[' mu_z = ', num2str(MODEL.mu_z),'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index b7d245783b9f132202b0141a76bddcdc2bfd6584..2c47363b5b599aced6933fd97b7f0da944b379b6 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 02;
+JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
data.FIGDIR = LOCALDIR;
@@ -38,24 +38,24 @@ if 0
options.T = [200 400];
fig = statistical_transport_averaging(data,options);
end
-
+
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
% options.NAME = '\phi';
-% options.NAME = 'N_i^{00}';
+options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
-options.NAME = '\Gamma_x';
-% options.NAME = 'n_i';
-options.PLAN = 'xy';
+% options.NAME = '\Gamma_x';
+% options.NAME = 'n_e';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 1;
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D(450:end);
+options.TIME = data.Ts3D(50:2:end);
% options.TIME = [350:600];
data.EPS = 0.1;
data.a = data.EPS * 2000;
@@ -69,16 +69,16 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = '\phi';
-% options.NAME = 'n_i';
-% options.NAME = 'N_e^{00}';
+options.NAME = 'n_e';
+% options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
-options.NAME = '\Gamma_x';
+% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [100 200 300 900 1300];
+options.TIME = [600];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 4fc0cb770121fb5ff2b21e9ab1bd5e78c778bfa2..815a21ade423539ba877a649f5edb8875a2324a5 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -10,13 +10,13 @@
% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
-% folder = '/misc/gene_results/ZP_16x8_kn_1.6/';
-folder = '/misc/gene_results/ZP_HP_kn_2.5/';
+folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
+% folder = '/misc/gene_results/ZP_HP_kn_2.5/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.2*gene_data.Ts3D(end);
+options.TAVG_0 = 0.6*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)
@@ -35,22 +35,22 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = 'Q_x';
options.NAME = '\phi';
-% options.NAME = 'T_i';
+% options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [0];
+options.TIME = [600];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
-save_figure(gene_data,fig)
+save_figure(gene_data,fig,'.png')
end
if 0
diff --git a/wk/g2k3_transport_scaling.m b/wk/g2k3_transport_scaling.m
index 8e83776846dc88aa3b62a0b13df4c080f1326fdc..3f06c1a2a52b62f696e30755b82d4362f1b858d0 100644
--- a/wk/g2k3_transport_scaling.m
+++ b/wk/g2k3_transport_scaling.m
@@ -29,7 +29,7 @@ figure;
% plot(K_N,G_LR); hold on;
% plot(K_N,G_LD); hold on;
% plot(K_N,G_CL,'--k'); hold on;
-k = 0.3;
+k = 0.35;
plot(K_N,g_SG.^2/k.^3); hold on;
plot(K_N,g_DG.^2/k.^3); hold on;
plot(K_N,g_LR.^2/k.^3); hold on;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index e703b67870e549b9a99522bed2649bb4f00b0c1c..3b3353860062c7a1a6631a065bed43af52cbe743 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -28,5 +28,8 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
%% ZPINCH
% outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3';
-outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
+outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
+% outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
+% outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
+% outfile ='Zpinch_rerun/Kn_1.6_200x48x21x11';
run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 55cc3c543eef6517bcec1734dea4115191587e42..687ef2dcc248b4d9e567265f3fe0c05408843e01 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -15,21 +15,21 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.0; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 2.5; % Density gradient drive
-K_T = 0.25*K_N; % Temperature '''
+K_N = 1.0;%2.0; % Density gradient drive
+K_T = 0;%0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
%% GRID PARAMETERS
-PMAXE = 4; % Hermite basis size of electrons
-JMAXE = 2; % Laguerre "
-PMAXI = 4; % " ions
-JMAXI = 2; % "
+PMAXE = 30; % Hermite basis size of electrons
+JMAXE = 15; % Laguerre "
+PMAXI = 30; % " ions
+JMAXI = 15; % "
NX = 1; % real space x-gridpoints
NY = 32; % '' y-gridpoints
LX = 100; % Size of the squared frequency domain
-LY = 60; % Size of the squared frequency domain
+LY = 30; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -41,7 +41,7 @@ SHEAR = 0.0; % magnetic shear (Not implemented yet)
EPS = 0.0; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 500; % Maximal time unit
-DT = 2e-2; % Time step
+DT = 5e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
@@ -53,7 +53,7 @@ SIMID = 'dbg'; % Name of the simulation
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'SG';
+CO = 'DG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -76,12 +76,13 @@ MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
+N_HD = 4;
MU_Z = 0.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
-NOISE0 = 0.0e-5; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
GRADB = 1.0;
CURVB = 1.0;
%%-------------------------------------------------------------------------