diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 5fa860f7b7b6e0c35479dd6f834a0a8dcbceaa59..b511eb92b17a8faa087253f5a88bcc910dd66746 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -42,7 +42,7 @@ while(CONTINUE)
W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
- W_NAPJ = 0;%strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
+ W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
W_SAPJ = 0;%strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
@@ -201,6 +201,7 @@ DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_;
DATA.PHI = PHI_;
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
+DATA.Je = Je; DATA.Ji = Ji;
DATA.kx = kx; DATA.ky = ky; DATA.z = z;
DATA.x = x; DATA.y = y;
DATA.ikx0 = ikx0; DATA.iky0 = iky0;
diff --git a/matlab/AssociatedLaguerrePoly.m b/matlab/compute/AssociatedLaguerrePoly.m
similarity index 100%
rename from matlab/AssociatedLaguerrePoly.m
rename to matlab/compute/AssociatedLaguerrePoly.m
diff --git a/matlab/compute/HermitePoly.m b/matlab/compute/HermitePoly.m
new file mode 100644
index 0000000000000000000000000000000000000000..e5c0bc5829c3dc8cf1977b5c12e6f8b3be122d26
--- /dev/null
+++ b/matlab/compute/HermitePoly.m
@@ -0,0 +1,40 @@
+
+% HermitePoly.m by David Terr, Raytheon, 5-10-04
+
+% Given nonnegative integer n, compute the
+% Hermite polynomial H_n. Return the result as a vector whose mth
+% element is the coefficient of x^(n+1-m).
+% polyval(HermitePoly(n),x) evaluates H_n(x).
+
+
+function hk = HermitePoly(n)
+% Evaluate the normalized Hermite polynomial.
+if n==0
+ hk = 1;
+elseif n==1
+ hk = [2 0];
+else
+
+ hkm2 = zeros(1,n+1);
+ hkm2(n+1) = 1;
+ hkm1 = zeros(1,n+1);
+ hkm1(n) = 2;
+
+ for k=2:n
+
+ hk = zeros(1,n+1);
+
+ for e=n-k+1:2:n
+ hk(e) = 2*(hkm1(e+1) - (k-1)*hkm2(e));
+ end
+
+ hk(n+1) = -2*(k-1)*hkm2(n+1);
+
+ if k<n
+ hkm2 = hkm1;
+ hkm1 = hk;
+ end
+
+ end
+
+end
\ No newline at end of file
diff --git a/matlab/LaguerrePoly.m b/matlab/compute/LaguerrePoly.m
similarity index 100%
rename from matlab/LaguerrePoly.m
rename to matlab/compute/LaguerrePoly.m
diff --git a/matlab/LinearFit_s.m b/matlab/compute/LinearFit_s.m
similarity index 100%
rename from matlab/LinearFit_s.m
rename to matlab/compute/LinearFit_s.m
diff --git a/matlab/computeFXYZ.m b/matlab/compute/computeFXYZ.m
similarity index 100%
rename from matlab/computeFXYZ.m
rename to matlab/compute/computeFXYZ.m
diff --git a/matlab/compute/compute_fa.m b/matlab/compute/compute_fa.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d8e4ec4ca83834a6d624f42738216601da4e375
--- /dev/null
+++ b/matlab/compute/compute_fa.m
@@ -0,0 +1,77 @@
+function [SS,XX,FF] = compute_fa(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));
+
+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);
+ otherwise
+ [~,iz] = min(abs(options.Z-data.z));
+ Napj_ = Napj_(:,:,:,:,iz,:);
+end
+Napj_ = squeeze(Napj_);
+
+frames = options.T;
+for it = 1:numel(options.T)
+ [~,frames(it)] = min(abs(options.T(it)-data.Ts5D));
+end
+
+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);
+ 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;
+ end
+ 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)
+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/compute/compute_fa_1D.m b/matlab/compute/compute_fa_1D.m
new file mode 100644
index 0000000000000000000000000000000000000000..3fce392234b36857c8aee44554bcf746313dd3f3
--- /dev/null
+++ b/matlab/compute/compute_fa_1D.m
@@ -0,0 +1,90 @@
+function [s,x, Fs, Fx] = compute_fa_1D(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);
+
+s = options.SPAR;
+x = options.XPERP;
+smin = min(abs(s));
+xmin = min(abs(x));
+
+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);
+ otherwise
+ [~,iz] = min(abs(options.Z-data.z));
+ Napj_ = Napj_(:,:,:,:,iz,:);
+end
+Napj_ = squeeze(Napj_);
+
+frames = options.T;
+for it = 1:numel(options.T)
+ [~,frames(it)] = min(abs(options.T(it)-data.Ts5D));
+end
+
+Napj_ = mean(Napj_(:,:,:,:,frames),5);
+
+Napj_ = squeeze(Napj_);
+
+
+Np = numel(parray); Nj = numel(jarray);
+
+% x = 0
+Fs = zeros(data.Nkx,data.Nky,numel(s));
+FAM = FaM(s,xmin);
+for ip_ = 1:Np
+ p_ = parray(ip_);
+ HH = Hp(p_,s);
+ for ij_ = 1:Nj
+ j_ = jarray(ij_);
+ LL = Lj(j_,xmin);
+ HLF = HH.*LL.*FAM;
+ for ikx = 1:data.Nkx
+ for iky = 1:data.Nky
+ Fs(ikx,iky,:) = squeeze(Fs(ikx,iky,:))' + Napj_(ip_,ij_,ikx,iky)*HLF;
+ end
+ end
+ end
+end
+
+% s = 0
+Fx = zeros(data.Nkx,data.Nky,numel(x));
+FAM = FaM(smin,x);
+for ip_ = 1:Np
+ p_ = parray(ip_);
+ HH = Hp(p_,smin);
+ for ij_ = 1:Nj
+ j_ = jarray(ij_);
+ LL = Lj(j_,x);
+ HLF = HH.*LL.*FAM;
+ for ikx = 1:data.Nkx
+ for iky = 1:data.Nky
+ Fx(ikx,iky,:) = squeeze(Fx(ikx,iky,:))' + Napj_(ip_,ij_,ikx,iky)*HLF;
+ end
+ end
+ end
+end
+
+Fs = real(Fs.*conj(Fs)); %|f_a|^2
+Fs = sqrt(squeeze(sum(sum(Fs,1),2))); %sqrt(<|f_a|>kx,ky)
+Fs = Fs./max(max(Fs));
+Fx = real(Fx.*conj(Fx)); %|f_a|^2
+Fx = sqrt(squeeze(sum(sum(Fx,1),2))); %sqrt(<|f_a|>kx,ky)
+Fx = Fx./max(max(Fx));
+end
\ No newline at end of file
diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
new file mode 100644
index 0000000000000000000000000000000000000000..01cf310596f9fb43e54fca71a6b80c3e218b805c
--- /dev/null
+++ b/matlab/compute/compute_fa_2D.m
@@ -0,0 +1,77 @@
+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));
+
+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);
+ otherwise
+ [~,iz] = min(abs(options.Z-data.z));
+ Napj_ = Napj_(:,:,:,:,iz,:);
+end
+Napj_ = squeeze(Napj_);
+
+frames = options.T;
+for it = 1:numel(options.T)
+ [~,frames(it)] = min(abs(options.T(it)-data.Ts5D));
+end
+
+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);
+ 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;
+ end
+ 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)
+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/dnjs.m b/matlab/compute/dnjs.m
similarity index 100%
rename from matlab/dnjs.m
rename to matlab/compute/dnjs.m
diff --git a/matlab/dnjs_fact.m b/matlab/compute/dnjs_fact.m
similarity index 100%
rename from matlab/dnjs_fact.m
rename to matlab/compute/dnjs_fact.m
diff --git a/matlab/dnjs_stir.m b/matlab/compute/dnjs_stir.m
similarity index 100%
rename from matlab/dnjs_stir.m
rename to matlab/compute/dnjs_stir.m
diff --git a/matlab/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m
similarity index 100%
rename from matlab/ifourier_GENE.m
rename to matlab/compute/ifourier_GENE.m
diff --git a/matlab/kernel.m b/matlab/compute/kernel.m
similarity index 100%
rename from matlab/kernel.m
rename to matlab/compute/kernel.m
diff --git a/matlab/compute_fa.m b/matlab/compute_fa.m
deleted file mode 100644
index 750cfc20f74db0a889ffef518ac296b778b08032..0000000000000000000000000000000000000000
--- a/matlab/compute_fa.m
+++ /dev/null
@@ -1,23 +0,0 @@
-function FF = compute_fa(Napj, spar, xperp)
-%% Compute the dispersion function from the moment hierarchi decomp.
-% Normalized Hermite
-Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p));
-% Laguerre
-Lj = @(j,x) polyval(LaguerrePoly(j),x);
-% Maxwellian factor
-FaM = @(s,x) exp(-s.^2-x);
-
-[SS, XX] = meshgrid(spar, xperp); %meshgrid for efficient evaluation
-
-FF = 0 .* SS;
-
-[Pmax,Jmax] = size(Napj);
-
-FAM = FaM(SS,XX);
-for p_ = 0:Pmax-1
- HH = Hp(p_,SS);
- for j_ = 0:Jmax-1
- LL = Lj(j_,XX);
- FF = FF + Napj(p_+1,j_+1)*HH.*LL.*FAM;
- end
-end
\ No newline at end of file
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 1dd71204e95f7d0c1a756d4778c733587667c2e7..e82fcf97e54f600d215cca10f1f93a3d8f5be810 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -1,10 +1,38 @@
fig = gcf;
axObjs = fig.Children;
dataObjs = axObjs.Children;
-
-X_ = dataObjs(1).XData;
-Y_ = dataObjs(1).YData;
-n0 = 0001;
+X_ = [];
+Y_ = [];
+for i = 1:numel(dataObjs)
+ X_ = [X_ dataObjs(i).XData];
+ Y_ = [Y_ dataObjs(i).YData];
+end
+n0 = 100;
figure;
-plot(X_(n0:end),Y_(n0:end));
-plot(X_(n0:end)-X_(n0),Y_(n0:end));
\ No newline at end of file
+mvm = @(x) movmean(x,1);
+shift = X_(n0);
+% shift = 0;
+% plot(X_(n0:end),Y_(n0:end));
+plot(mvm(X_(n0:end)),mvm(Y_(n0:end))./max(Y_(200:end))); hold on;
+
+
+n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
+avg_ = mean(Y_(n1:n2));
+std_ = std(Y_(n1:n2));
+title(['avg = ',num2str(avg_),' std = ', num2str(std_)])
+plot([X_(n1),X_(n2)],[1 1]*avg_,'--k')
+
+% ylim([0,avg_*3]);
+sum(Y_(2:end)./X_(2:end).^(3/2))
+%%
+
+if 0
+ %%
+ m1 = mean(y1); m2 = mean(y2);
+ s1 = std(y1); s2 = std(y2);
+ delta = abs(m2-m1)/min(s1,s2);
+ delta
+end
+
+
+ xlim([-3,3]); ylim([0,4.5]);
\ No newline at end of file
diff --git a/matlab/loadGeneral_GENE.m b/matlab/load/loadGeneral_GENE.m
similarity index 100%
rename from matlab/loadGeneral_GENE.m
rename to matlab/load/loadGeneral_GENE.m
diff --git a/matlab/load_0D_data.m b/matlab/load/load_0D_data.m
similarity index 100%
rename from matlab/load_0D_data.m
rename to matlab/load/load_0D_data.m
diff --git a/matlab/load_2D_data.m b/matlab/load/load_2D_data.m
similarity index 100%
rename from matlab/load_2D_data.m
rename to matlab/load/load_2D_data.m
diff --git a/matlab/load_3D_data.m b/matlab/load/load_3D_data.m
similarity index 100%
rename from matlab/load_3D_data.m
rename to matlab/load/load_3D_data.m
diff --git a/matlab/load_4D_data.m b/matlab/load/load_4D_data.m
similarity index 100%
rename from matlab/load_4D_data.m
rename to matlab/load/load_4D_data.m
diff --git a/matlab/load_5D_data.m b/matlab/load/load_5D_data.m
similarity index 100%
rename from matlab/load_5D_data.m
rename to matlab/load/load_5D_data.m
diff --git a/matlab/load_daint.m b/matlab/load/load_daint.m
similarity index 100%
rename from matlab/load_daint.m
rename to matlab/load/load_daint.m
diff --git a/matlab/load_field_gene.m b/matlab/load/load_field_gene.m
similarity index 71%
rename from matlab/load_field_gene.m
rename to matlab/load/load_field_gene.m
index a323a9a277240a7dc0a28dd5f6af36eff4b0256f..8095418eadce4de729d0c761d52aa6135ce1b550 100644
--- a/matlab/load_field_gene.m
+++ b/matlab/load/load_field_gene.m
@@ -1,7 +1,10 @@
% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_KHR_2/';
% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-3_gyroLES/';
% folder = '/misc/gene_results/HP_fig_2c_gyroLES/';
-folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK/';
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK/';
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_30x30_Lv_smaller/';
+% folder = '/misc/gene_results/HP_fig_2b_mu_1e-1/';
+folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/';
%% General info (domain etc.)
varargin{1} = 0; varargin{2} = 0;
@@ -17,7 +20,7 @@ x = Lx/nx.*(-nx/2:nx/2-1);
y = Ly/ny.*(-ny+1:ny-1);
%% Load dens_i
-steps_range = 0:1:500; nt = numel(steps_range);
+steps_range = 0:1:200; nt = numel(steps_range);
time_dens = 1:nt;
dens_xt = zeros(nt,nx);
% figure
@@ -33,7 +36,7 @@ end
if 1
%% Load phi
-steps_range = 0:2:1000; nt = numel(steps_range);
+steps_range = 0:1:1000; nt = numel(steps_range);
time_phi = 1:nt;
phi_xt = zeros(nt,nx);
% figure
@@ -64,7 +67,8 @@ for i = 1:nt
time_vExB(i) = t;
end
end
-if 1
+
+if 0
%% compute Gx
steps_n = 0:1:500; nt1 = numel(steps_n);
steps_v = 0:2:1000; nt2 = numel(steps_v);
@@ -85,14 +89,36 @@ for i = 1:nt1
Gx(i) = mean(mean(squeeze(dens(:,:,1).*vExB(:,:,1))));
end
end
-%% Plot
+%% Plot spacetime
figure;
-[X_TX,Y_TX] = meshgrid(time_dens,x);
+% [X_TX,Y_TX] = meshgrid(time_dens,x);
% pclr=pcolor(X_TX',Y_TX',dens_xt);
+[X_TX,Y_TX] = meshgrid(time_phi,x);
pclr=pcolor(X_TX',Y_TX',phi_xt);
set(pclr, 'edgecolor','none');
colormap(bluewhitered(256))
-%%
-% figure
-% plot(time_dens,Gx)
+%% Load <f>
+
+fid=fopen([folder,'vsp.dat']);
+
+nz=general.allParams.box.nz;
+nv=general.allParams.box.nv;
+nw=general.allParams.box.nw;
+ns=general.allParams.specs.ns;
+
+[~,~,gene_vsp]=GetFortranStep(fid,'DOUBLE',8,'LITTLE',1);
+fclose(fid);
+
+ gridsize=prod([nz nv nw]);
+
+
+%this is the least effective way of doing it, but fine
+i_s = 1;
+shift=(i_s-1)*gridsize;
+f_avg=reshape(gene_vsp(shift+4*gridsize+1:shift+5*gridsize),[nz nv nw]);
+
+figure
+pclr=pcolor(squeeze(f_avg(1,:,:))')
+set(pclr, 'edgecolor','none');
+colormap(bluewhitered(256))
diff --git a/matlab/load_grid_data.m b/matlab/load/load_grid_data.m
similarity index 100%
rename from matlab/load_grid_data.m
rename to matlab/load/load_grid_data.m
diff --git a/matlab/load_marconi.m b/matlab/load/load_marconi.m
similarity index 100%
rename from matlab/load_marconi.m
rename to matlab/load/load_marconi.m
diff --git a/matlab/load_multiple_outputs_marconi.m b/matlab/load/load_multiple_outputs_marconi.m
similarity index 100%
rename from matlab/load_multiple_outputs_marconi.m
rename to matlab/load/load_multiple_outputs_marconi.m
diff --git a/matlab/load_params.m b/matlab/load/load_params.m
similarity index 100%
rename from matlab/load_params.m
rename to matlab/load/load_params.m
diff --git a/matlab/load_results.m b/matlab/load/load_results.m
similarity index 100%
rename from matlab/load_results.m
rename to matlab/load/load_results.m
diff --git a/matlab/plot/1D_time_average_spectrum.m b/matlab/plot/1D_time_average_spectrum.m
new file mode 100644
index 0000000000000000000000000000000000000000..d776e9ce42061e83a179ecf8d9e57b1dda364687
--- /dev/null
+++ b/matlab/plot/1D_time_average_spectrum.m
@@ -0,0 +1,69 @@
+function [ FIGURE ] = spectrum_1D( data, options )
+
+options.PLAN = 'kxky';
+options.COMP = 1;
+
+toplot = process_field(data,options);
+t = data.Ts3D; frames = toplot.FRAMES;
+
+colors = jet(numel(frames));
+
+switch options.NAME
+ case '\Gamma_x';
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS];
+ case '\phi'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS];
+ case 'n_i'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS];
+ case 'n_e'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS];
+end
+set(gcf, 'Position', [20 50 5000 2000])
+subplot(1,2,1)
+
+ sdim = 2;
+ k = data.kx;
+ xname = '$k_x$';
+ yname = '$\sum_{k_y}|\Gamma_k|$';
+ nmax = ceil(data.Nkx*2/3);
+ shiftx = @(x) x(1:nmax);
+ shifty = @(x) x(1:nmax);
+
+ for it = 1:numel(toplot.FRAMES)
+ Gk = sum(abs(toplot.FIELD(:,:,toplot.FRAMES(it))),sdim);
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ semilogy(shiftx(k), shifty(Gk),'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ grid on
+ title('HeLaZ $k_x$ transport spectrum'); legend('show','Location','eastoutside')
+ xlabel(xname); ylabel(yname)
+
+subplot(1,2,2)
+
+ sdim = 1;
+ k = data.ky;
+ xname = '$k_y$';
+ yname = '$\sum_{k_x}|\Gamma_k|$';
+ nmax = floor(data.Nky/2*2/3);
+ AA = @(x) x(1:nmax);
+ shiftx = @(x) AA(x);
+ shifty = @(x) AA(ifftshift(x));
+
+ for it = 1:numel(toplot.FRAMES)
+ Gk = sum(abs(toplot.FIELD(:,:,toplot.FRAMES(it))),sdim);
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ semilogy(shiftx(k), shifty(Gk),'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ grid on
+ title('HeLaZ $k_y$ transport spectrum'); legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(yname)
+end
+
diff --git a/matlab/ZF_fourier_analysis.m b/matlab/plot/ZF_fourier_analysis.m
similarity index 100%
rename from matlab/ZF_fourier_analysis.m
rename to matlab/plot/ZF_fourier_analysis.m
diff --git a/matlab/plots/ZF_spacetime.m b/matlab/plot/ZF_spacetime.m
similarity index 100%
rename from matlab/plots/ZF_spacetime.m
rename to matlab/plot/ZF_spacetime.m
diff --git a/matlab/bluewhitered.m b/matlab/plot/bluewhitered.m
similarity index 100%
rename from matlab/bluewhitered.m
rename to matlab/plot/bluewhitered.m
diff --git a/matlab/create_film.m b/matlab/plot/create_film.m
similarity index 64%
rename from matlab/create_film.m
rename to matlab/plot/create_film.m
index 5f9724726d2821d1ba205ca01925bc59a4e6d225..ed66da62b4d2d4ea61c226f29b25bda5870913cb 100644
--- a/matlab/create_film.m
+++ b/matlab/plot/create_film.m
@@ -2,7 +2,11 @@ function create_film(DATA,OPTIONS,format)
%% Plot options
FPS = 30; DELAY = 1/FPS;
BWR = 1; NORMALIZED = 1;
-T = DATA.Ts3D;
+if ~strcmp(OPTIONS.PLAN,'sx')
+ T = DATA.Ts3D;
+else
+ T = DATA.Ts5D;
+end
%% Processing
toplot = process_field(DATA,OPTIONS);
@@ -29,40 +33,49 @@ if hmax == hmin
else
% Setup figure frame
fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]);
- pcolor(X,Y,FIELD(:,:,1)); % to set up
- if BWR
- colormap(bluewhitered)
- end
- axis tight manual % this ensures that getframe() returns a consistent size
- if toplot.INTERP
- shading interp;
+ if ~strcmp(OPTIONS.PLAN,'sx')
+ pcolor(X,Y,FIELD(:,:,1)); % to set up
+ if BWR
+ colormap(bluewhitered)
+ end
+ axis tight manual % this ensures that getframe() returns a consistent size
+ if toplot.INTERP
+ shading interp;
+ end
+ else
+ contour(toplot.X,toplot.Y,FIELD(:,:,1),128);
end
in = 1;
nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
for n = FRAMES % loop over selected frames
scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize
- if NORMALIZED
- pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
- caxis([-1,1]);
- else
- pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
- if CONST_CMAP
- caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map
+ if ~strcmp(OPTIONS.PLAN,'sx')
+ if NORMALIZED
+ pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
+ caxis([-1,1]);
else
- caxis([-1,1]*scale); % adaptive color map
+ pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
+ if CONST_CMAP
+ caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map
+ else
+ caxis([-1,1]*scale); % adaptive color map
+ end
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scale = ',sprintf('%.1e',scale)]);
end
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scale = ',sprintf('%.1e',scale)]);
+ if toplot.INTERP
+ shading interp;
+ end
+ set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT);
+ if BWR
+ colormap(bluewhitered)
+ end
+ else % show velocity distr.
+ contour(toplot.X,toplot.Y,FIELD(:,:,n)/scale,128);
end
title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
,', scaling = ',sprintf('%.1e',scale)]);
- if toplot.INTERP
- shading interp;
- end
- set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT);
- if BWR
- colormap(bluewhitered)
- end
+
xlabel(XNAME); ylabel(YNAME); %colorbar;
drawnow
% Capture the plot as an image
diff --git a/matlab/default_plots_options.m b/matlab/plot/default_plots_options.m
similarity index 100%
rename from matlab/default_plots_options.m
rename to matlab/plot/default_plots_options.m
diff --git a/matlab/distinguishable_colors.m b/matlab/plot/distinguishable_colors.m
similarity index 100%
rename from matlab/distinguishable_colors.m
rename to matlab/plot/distinguishable_colors.m
diff --git a/matlab/plot/gene_transport_spectrum.m b/matlab/plot/gene_transport_spectrum.m
new file mode 100644
index 0000000000000000000000000000000000000000..66563af54d6fb57481d20b6b195b3ec53f2a03c4
--- /dev/null
+++ b/matlab/plot/gene_transport_spectrum.m
@@ -0,0 +1,99 @@
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_36x20/';
+% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
+folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/HP_fig_2c_gyroLES/';
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_128x32x48x32/';
+
+file = 'coord.dat.h5';
+vp = h5read([folder,file],'/coord/vp');
+mu = h5read([folder,file],'/coord/mu');
+kx = h5read([folder,file],'/coord/kx');
+ky = h5read([folder,file],'/coord/ky');
+z = h5read([folder,file],'/coord/z');
+
+[KY,~] = meshgrid(ky,kx);
+
+file = 'field.dat.h5';
+time = h5read([folder,file],'/field/time');
+
+KDIM = 'ky';
+TIMES = 2500:100:3500;
+switch KDIM
+ case 'kx'
+ sdim = 2;
+ k = kx;
+ xname = '$k_x$';
+ yname = '$\sum_{k_y}|\Gamma_k|$';
+ shiftx = @(x) x(1:numel(kx)/2);
+ shifty = @(x) x(1:numel(kx)/2);
+ case 'ky'
+ sdim = 1;
+ k = ky;
+ xname = '$k_y$';
+ yname = '$\sum_{k_x}|\Gamma_k|$';
+ shiftx = @(x) x;
+ shifty = @(x) x(1:numel(ky));
+end
+
+colors = jet(numel(TIMES));
+i_ = 1;
+figure
+for T = TIMES
+[~, it] = min(abs(time-T));
+
+v_x = h5read([folder, 'field.dat.h5'],[ '/field/phi/',sprintf('%10.10d',it-1)]);
+v_x = v_x.real + 1i*v_x.imaginary; v_x = squeeze(v_x(:,:,1));
+v_x = -1i*KY.*v_x;
+
+n_i = h5read([folder,'mom_ions.dat.h5'],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+n_i = n_i.real + 1i*n_i.imaginary; n_i = squeeze(n_i(:,:,1));
+% topclr_ = abs(n_i);
+
+% v_x = computeFXYZ(v_x);
+% n_i = computeFXYZ(n_i);
+v_x = real(ifft2(v_x));
+n_i = real(ifft2(n_i));
+
+Gx = n_i.*v_x;
+
+Gx = fft2(Gx);
+% topclr_ = abs(Gx);
+end
+[~, it0] = min(abs(time-T(1)));
+[~, it1] = min(abs(time-T(end)));
+Gx = mean(Gx,3);
+
+subplot(1,2,1)
+ sdim = 2;
+ k = kx;
+ shiftx = @(x) x;%(1:numel(kx)/2);
+ shifty = @(x) x;%(1:numel(kx)/2);
+ Y = abs(Gx(:,ky==0));
+ Y = Y/max(Y);
+ semilogy(shiftx(k),shifty(Y),'DisplayName',['$t=',num2str(time(it)),'$'],...
+ 'Color',colors(i_,:)); hold on;
+
+subplot(1,2,2)
+ sdim = 1;
+ k = ky;
+ shiftx = @(x) x;
+ shifty = @(x) x;%(1:numel(ky));
+ Y = abs(Gx(kx==0,:));
+ Y = Y/max(Y);
+ semilogy(shiftx(k),shifty(Y),'DisplayName',['$t=',num2str(time(it)),'$'],...
+ 'Color',colors(i_,:)); hold on;
+i_ = i_+1;
+
+subplot(1,2,1)
+ xname = '$k_x$';
+ yname = '$\sum_{k_y}|\Gamma_k|$';
+ title('Gene $k_x$ transport spectrum'); legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(yname)
+
+subplot(1,2,2)
+ xname = '$k_y$';
+ yname = '$\sum_{k_x}|\Gamma_k|$';
+ title('Gene $k_y$ transport spectrum'); legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(yname)
+
+ figure; pclr = pcolor(topclr_); set(pclr,'EdgeColor','none');
diff --git a/matlab/open_figure_script.m b/matlab/plot/open_figure_script.m
similarity index 100%
rename from matlab/open_figure_script.m
rename to matlab/plot/open_figure_script.m
diff --git a/matlab/photomaton.m b/matlab/plot/photomaton.m
similarity index 50%
rename from matlab/photomaton.m
rename to matlab/plot/photomaton.m
index 66b6f04832b2f43315fb7fe50570158fc60eb802..4a780f350a12ff72f87002d7e9157f4dcf51f6bb 100644
--- a/matlab/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -17,20 +17,28 @@ FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
else
scale = 1;
end
- pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none');
- if OPTIONS.AXISEQUAL
- pbaspect(toplot.ASPECT)
- end
- if ~strcmp(OPTIONS.PLAN,'kxky')
- caxis([-1,1]);
- colormap(bluewhitered);
+ if ~strcmp(OPTIONS.PLAN,'sx')
+ tshot = DATA.Ts3D(FRAMES(i_));
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none');
+ if OPTIONS.AXISEQUAL
+ pbaspect(toplot.ASPECT)
+ end
+ if ~strcmp(OPTIONS.PLAN,'kxky')
+ caxis([-1,1]);
+ colormap(bluewhitered);
+ end
+ if OPTIONS.INTERP
+ shading interp;
+ end
+ else
+ contour(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale,128);
+% pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none'); shading interp
+ tshot = DATA.Ts5D(FRAMES(i_));
end
xlabel(toplot.XNAME); ylabel(toplot.YNAME);
% if i_ > 1; set(gca,'ytick',[]); end;
- title([sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_))),', max = ',sprintf('%.1e',scale)]);
- if OPTIONS.INTERP
- shading interp;
- end
+ title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',scale)]);
+
end
legend(['$',toplot.FIELDNAME,'$']);
FIGURE.FIGNAME = [FNAME,'_snaps',TNAME];
diff --git a/matlab/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
similarity index 100%
rename from matlab/plot_cosol_mat.m
rename to matlab/plot/plot_cosol_mat.m
diff --git a/matlab/plot/plot_fa.m b/matlab/plot/plot_fa.m
new file mode 100644
index 0000000000000000000000000000000000000000..d95e2c44cdd8793863724b13d20ea48d1b6af00f
--- /dev/null
+++ b/matlab/plot/plot_fa.m
@@ -0,0 +1,60 @@
+function [ FIGURE ] = plot_fa( DATA, OPTIONS )
+
+FIGURE.fig = figure; FIGURE.FIGNAME = ['f_a_',DATA.PARAMS];
+
+if OPTIONS.ONED
+ OPTIONS.SPECIE = 'i';
+ [~,~,fsi,fxi] = compute_fa_1D(DATA, OPTIONS);
+ OPTIONS.SPECIE = 'e';
+ [s,x,fse,fxe] = compute_fa_1D(DATA, OPTIONS);
+ [~,it] = min(abs(OPTIONS.T-DATA.Ts5D));
+ subplot(1,2,1)
+ plot(s,fse); hold on
+ plot(s,fsi);
+ legend('e','i')
+ xlabel('$v_\parallel, (\mu=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+ title(DATA.param_title);
+ subplot(1,2,2)
+ plot(x,fxe); hold on;
+ plot(x,fxi);
+ legend('e','i')
+ xlabel('$\mu, (v_\parallel=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+ if numel(it) == 1
+ title(['t=',num2str(DATA.Ts5D(it))]);
+ else
+ title(['average $t\in$[',num2str(DATA.Ts5D(min(it))),',',num2str(DATA.Ts5D(max(it))),']']);
+ end
+
+else
+ OPTIONS.SPECIE = 'i';
+ [~,~,FFi] = compute_fa_2D(DATA, OPTIONS);
+ OPTIONS.SPECIE = 'e';
+ [SS,XX,FFe] = compute_fa_2D(DATA, OPTIONS);
+ [~,it] = min(abs(OPTIONS.T-DATA.Ts5D));
+ subplot(1,2,1)
+ if OPTIONS.CTR
+ contour(SS,XX,FFi',128);
+ else
+ pclr = pcolor(SS,XX,FFi'); set(pclr, 'edgecolor','none'); shading interp
+ end
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
+ legend('$\langle |f_i|^2\rangle_{xy}^{1/2}$')
+ title(DATA.param_title);
+ % FF = log10(FF);
+ subplot(1,2,2)
+ if OPTIONS.CTR
+ contour(SS,XX,FFe',128);
+ else
+ pclr = pcolor(SS,XX,FFe'); set(pclr, 'edgecolor','none'); shading interp
+ end
+ legend('$\langle |f_e|^2\rangle_{xy}^{1/2}$')
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
+ if numel(it) == 1
+ title(['t=',num2str(DATA.Ts5D(it))]);
+ else
+ title(['average $t\in$[',num2str(DATA.Ts5D(min(it))),',',num2str(DATA.Ts5D(max(it))),']']);
+ end
+end
+
+end
+
diff --git a/matlab/plot_kernels.m b/matlab/plot/plot_kernels.m
similarity index 100%
rename from matlab/plot_kernels.m
rename to matlab/plot/plot_kernels.m
diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plot/plot_kperp_spectrum.m
similarity index 100%
rename from matlab/plots/plot_kperp_spectrum.m
rename to matlab/plot/plot_kperp_spectrum.m
diff --git a/matlab/plot_param_evol.m b/matlab/plot/plot_param_evol.m
similarity index 64%
rename from matlab/plot_param_evol.m
rename to matlab/plot/plot_param_evol.m
index 827b4c931985bf1db60ba83071f215b392152840..b8707232eacfe13883a0a618240edb278fc1d6ef 100644
--- a/matlab/plot_param_evol.m
+++ b/matlab/plot/plot_param_evol.m
@@ -1,16 +1,28 @@
-fig = figure; FIGNAME = 'linear_study';
+function [fig] = plot_param_evol(data)
+TJOB_SE = data.TJOB_SE;
+NU_EVOL = data.NU_EVOL;
+DT_EVOL = data.DT_EVOL;
+CO_EVOL = data.CO_EVOL;
+MUx_EVOL = data.MUx_EVOL;
+MUy_EVOL = data.MUy_EVOL;
+% K_T_EVOL = data.K_T_EVOL;
+K_N_EVOL = data.K_N_EVOL;
+L_EVOL = data.L_EVOL;
+
+fig = figure;
hold on; set(gcf, 'Position', [100, 100, 400, 1500])
subplot(4,1,1);
title('Parameter evolution'); hold on;
yyaxis left
plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
yyaxis right
- plot(TJOB_SE,CO_EVOL,'DisplayName','CO');
+ plot(TJOB_SE,DT_EVOL,'DisplayName','dt');
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
plot_tjob_lines(TJOB_SE,ylim)
subplot(4,1,2);
- plot(TJOB_SE,MU_EVOL,'DisplayName','$\mu$'); hold on;
+ plot(TJOB_SE,MUx_EVOL,'DisplayName','$\mu_x$'); hold on;
+ plot(TJOB_SE,MUy_EVOL,'DisplayName','$\mu_y$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
plot_tjob_lines(TJOB_SE,ylim)
@@ -19,7 +31,7 @@ subplot(4,1,3);
yyaxis left
plot(TJOB_SE,K_N_EVOL,'DisplayName','$\nabla n$'); hold on;
yyaxis right
- plot(TJOB_SE,K_T_EVOL,'DisplayName','$\nabla T$'); hold on;
+% plot(TJOB_SE,K_T_EVOL,'DisplayName','$\nabla T$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
plot_tjob_lines(TJOB_SE,ylim)
@@ -30,10 +42,11 @@ subplot(4,1,4);
plot_tjob_lines(TJOB_SE,ylim)
xlabel('$t c_s/R$');
- saveas(fig,[BASIC.RESDIR,'param_evol.png']);
function [] = plot_tjob_lines(TJOB,limits)
for i = 2:numel(TJOB)-1
plot(TJOB(i)*[1 1],limits,'--k')
end
end
+
+end
\ No newline at end of file
diff --git a/matlab/plot_phi_ballooning.m b/matlab/plot/plot_phi_ballooning.m
similarity index 100%
rename from matlab/plot_phi_ballooning.m
rename to matlab/plot/plot_phi_ballooning.m
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plot/plot_radial_transport_and_shear.m
similarity index 100%
rename from matlab/plots/plot_radial_transport_and_shear.m
rename to matlab/plot/plot_radial_transport_and_shear.m
diff --git a/matlab/plots/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
similarity index 87%
rename from matlab/plots/plot_radial_transport_and_spacetime.m
rename to matlab/plot/plot_radial_transport_and_spacetime.m
index e6039a696796b110a20f78221903092d9c1c1288..fa88f634363b26a4fb3f0288a339c46eebf39693 100644
--- a/matlab/plots/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,4 +1,4 @@
-function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname)
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname,Nmvm)
%Compute steady radial transport
tend = TAVG_1; tstart = TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
@@ -34,16 +34,17 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
[~,ite3D] = min(abs(DATA.Ts3D-tend));
%% Figure
+mvm = @(x) movmean(x,Nmvm);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
% yyaxis left
- plot(DATA.Ts0D,DATA.PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
- plot(DATA.Ts3D,Gamma_t_mtlb,'DisplayName','matlab comp.'); hold on;
+ plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+ plot(mvm(DATA.Ts3D),mvm(Gamma_t_mtlb),'DisplayName','matlab comp.'); hold on;
plot(DATA.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_x$')
ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
- title([DATA.param_title,', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
+ title([DATA.param_title,', $\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
% plot
subplot(312)
it0 = 1; itend = Ns3D;
@@ -64,6 +65,9 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
Ns3D = numel(DATA.Ts3D);
[KY, KX] = meshgrid(DATA.ky, DATA.kx);
plt = @(x) mean(x(:,:,1,:),2);
+ kycut = max(DATA.ky);
+ kxcut = max(DATA.kx);
+ LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
switch stfname
case 'phi'
phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
@@ -74,19 +78,19 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
case 'v_y'
dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
end
f2plot = dxphi; fname = '$\langle \partial_x\phi\rangle_y$';
case 'v_x'
dyphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dyphi(:,:,1,it) = ifourier_GENE(1i*KY.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ dyphi(:,:,1,it) = ifourier_GENE(1i*KY.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
end
f2plot = dyphi; fname = '$\langle \partial_y\phi\rangle_y$';
case 'szf'
dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
end
f2plot = dx2phi; fname = '$\langle \partial_x^2\phi\rangle_y$';
end
diff --git a/matlab/plots/plot_shear_and_reynold_stress.m b/matlab/plot/plot_shear_and_reynold_stress.m
similarity index 100%
rename from matlab/plots/plot_shear_and_reynold_stress.m
rename to matlab/plot/plot_shear_and_reynold_stress.m
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plot/plot_space_time_diagrams.m
similarity index 100%
rename from matlab/plots/plot_space_time_diagrams.m
rename to matlab/plot/plot_space_time_diagrams.m
diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plot/plot_time_evolution_and_gr.m
similarity index 100%
rename from matlab/plots/plot_time_evolution_and_gr.m
rename to matlab/plot/plot_time_evolution_and_gr.m
diff --git a/matlab/plot/real_plot_1D.m b/matlab/plot/real_plot_1D.m
new file mode 100644
index 0000000000000000000000000000000000000000..5b7eda421d3a8c4941c3c952d5347a03cc430c09
--- /dev/null
+++ b/matlab/plot/real_plot_1D.m
@@ -0,0 +1,113 @@
+function [ FIGURE ] = real_plot_1D( data, options )
+
+options.PLAN = 'xy';
+options.COMP = options.COMPZ;
+options.POLARPLOT = 0;
+options.INTERP = 0;
+
+toplot = process_field(data,options);
+t = data.Ts3D; frames = toplot.FRAMES;
+
+colors = jet(numel(frames));
+
+switch options.NAME
+ case '\Gamma_x'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['transport_1D_avg_',data.PARAMS];
+ yname = '\Gamma';
+ case 'v_y'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_1D_avg_',data.PARAMS];
+ yname = 'v_y';
+ case 's_y'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZS_1D_avg_',data.PARAMS];
+ yname = 's_y';
+ case '\phi'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_1D_avg_',data.PARAMS];
+ yname = '\phi';
+ case 'n_i'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_1D_avg_',data.PARAMS];
+ yname = 'n_i';
+ case 'n_e'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_1D_avg_',data.PARAMS];
+ yname = 'n_e';
+end
+
+switch options.COMPXY
+ case 'avg'
+ compx = @(x) mean(x,1);
+ compy = @(x) mean(x,2);
+ ynamex = ['$\langle ',yname,'\rangle_x$'];
+ ynamey = ['$\langle ',yname,'\rangle_y$'];
+ otherwise
+ compx = @(x) x(1,:);
+ compy = @(x) x(:,1);
+ ynamex = ['$',yname,'(x=0)$'];
+ ynamey = ['$',yname,'(y=0)$'];
+end
+
+
+set(gcf, 'Position', [20 50 5000 2000])
+subplot(1,2,1)
+
+ X = data.x;
+ xname = '$x$';
+ switch options.COMPT
+ case 'avg'
+ Y = compy(mean(toplot.FIELD(:,:,:),3));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t\in[',num2str(t(frames(1))),',',num2str(t(frames(end))),']$'],...
+ 'Color',colors(1,:)); hold on;
+ legend('show')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Y = compy(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ legend('show','Location','eastoutside')
+ end
+ grid on
+ xlabel(xname); ylabel(ynamey)
+ xlim([min(X),max(X)]);
+
+subplot(1,2,2)
+
+ X = data.y;
+ xname = '$y$';
+ switch options.COMPT
+ case 'avg'
+ Y = compx(mean(toplot.FIELD(:,:,:),3));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t\in[',num2str(t(frames(1))),',',num2str(t(frames(end))),']$'],...
+ 'Color',colors(1,:)); hold on;
+ legend('show')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Y = compx(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ legend('show','Location','eastoutside')
+ end
+
+ grid on
+% title('HeLaZ $k_y$ transport spectrum');
+ legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(ynamex)
+ xlim([min(X),max(X)]);
+
+end
+
diff --git a/matlab/save_figure.m b/matlab/plot/save_figure.m
similarity index 100%
rename from matlab/save_figure.m
rename to matlab/plot/save_figure.m
diff --git a/matlab/show_geometry.m b/matlab/plot/show_geometry.m
similarity index 100%
rename from matlab/show_geometry.m
rename to matlab/plot/show_geometry.m
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
new file mode 100644
index 0000000000000000000000000000000000000000..5dc78faef956803fac446f07a073d504b44b5e6f
--- /dev/null
+++ b/matlab/plot/show_moments_spectrum.m
@@ -0,0 +1,142 @@
+function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
+
+
+Pi = DATA.Pi;
+Ji = DATA.Ji;
+Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
+Nipj = squeeze(Nipj);
+
+Pe = DATA.Pe;
+Je = DATA.Je;
+Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
+Nepj = squeeze(Nepj);
+
+FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS];
+set(gcf, 'Position', [100 50 1000 400])
+
+if ~OPTIONS.ST
+ t0 = OPTIONS.TIME(1);
+ t1 = OPTIONS.TIME(end);
+ [~,it0] = min(abs(t0-DATA.Ts5D));
+ [~,it1] = min(abs(t1-DATA.Ts5D));
+ Nipj = mean(Nipj(:,:,it0:it1),3);
+ Nepj = mean(Nepj(:,:,it0:it1),3);
+ if numel(OPTIONS.TIME) == 1
+ TITLE=['$t=',num2str(OPTIONS.TIME),'$'];
+ else
+ TITLE=['$t\in[',num2str(t0),',',num2str(t1),']$'];
+ end
+ Nipj = squeeze(Nipj);
+ Nepj = squeeze(Nepj);
+
+ ymini = min(Nipj); ymaxi = max(Nipj);
+ ymine = min(Nepj); ymaxe = max(Nepj);
+ ymin = min([ymini ymine]);
+ ymax = max([ymaxi ymaxe]);
+
+
+ subplot(121)
+ if ~OPTIONS.P2J
+ for ij = 1:numel(Ji)
+ name = ['$j=',num2str(Ji(ij)),'$'];
+ semilogy(Pi,Nipj(:,ij),'o-','DisplayName',name); hold on;
+ end
+ xlabel('$p$');
+ else
+ for ij = 1:numel(Ji)
+ name = ['$j=',num2str(Ji(ij)),'$'];
+ semilogy(Pi+2*Ji(ij),Nipj(:,ij),'o-','DisplayName',name); hold on;
+ end
+ xlabel('$p+2j$')
+ end
+ ylabel(['$\sum_{kx,ky}|N_i^{pj}|$']);
+ legend('show');
+ title([TITLE,' He-La ion spectrum']);
+
+ subplot(122)
+ if ~OPTIONS.P2J
+ for ij = 1:numel(Je)
+ name = ['$j=',num2str(Je(ij)),'$'];
+ semilogy(Pe,Nepj(:,ij),'o-','DisplayName',name); hold on;
+ end
+ xlabel('p');
+ else
+ for ij = 1:numel(Je)
+ name = ['$j=',num2str(Je(ij)),'$'];
+ semilogy(Pe+2*Je(ij),Nepj(:,ij),'o-','DisplayName',name); hold on;
+ end
+ xlabel('$p+2j$')
+ end
+ ylabel(['$\sum_{kx,ky}|N_e^{pj}|$']);
+ legend('show');
+ title([TITLE,' He-La elec. spectrum']);
+else
+ [JJ,PP] = meshgrid(Ji,Pi);
+ P2Ji = PP + 2*JJ;
+ % form an axis of velocity ordered moments
+ p2ji = unique(P2Ji);
+ % weights to average
+ weights = zeros(size(p2ji));
+ % space time of moments amplitude wrt to p+2j degree
+ Ni_ST = zeros(numel(p2ji),numel(DATA.Ts5D));
+ % fill the st diagramm by averaging every p+2j=n moments together
+ for ip = 1:numel(Pi)
+ for ij = 1:numel(Ji)
+ [~,ip2j] = min(abs(p2ji-(Pi(ip)+2*Ji(ij))));
+ Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
+ weights(ip2j) = weights(ip2j) + 1;
+ end
+ end
+ % doing the average
+ for ip2j = 1:numel(p2ji)
+ Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
+ end
+ % same for electrons!!
+ [JJ,PP] = meshgrid(Je,Pe);
+ P2Je = PP + 2*JJ;
+ % form an axis of velocity ordered moments
+ p2je = unique(P2Je);
+ % weights to average
+ weights = zeros(size(p2je));
+ % space time of moments amplitude wrt to p+2j degree
+ Ne_ST = zeros(numel(p2je),numel(DATA.Ts5D));
+ % fill the st diagramm by averaging every p+2j=n moments together
+ for ip = 1:numel(Pe)
+ for ij = 1:numel(Je)
+ [~,ip2j] = min(abs(p2ji-(Pe(ip)+2*Je(ij))));
+ Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
+ weights(ip2j) = weights(ip2j) + 1;
+ end
+ end
+ % doing the average
+ for ip2j = 1:numel(p2ji)
+ Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
+ end
+ % plots
+ % scaling
+% plt = @(x,ip2j) x./max(max(x));
+ if OPTIONS.NORMALIZED
+ plt = @(x,ip2j) x(ip2j,:)./max(x(ip2j,:));
+ else
+ plt = @(x,ip2j) x;
+ end
+ subplot(2,1,1)
+ imagesc(DATA.Ts5D,p2ji,plt(Ni_ST,1:numel(p2ji)));
+ set(gca,'YDir','normal')
+% pclr = pcolor(XX,YY,plt(Ni_ST));
+% set(pclr, 'edgecolor','none'); hold on;
+ xlabel('$t$'); ylabel('$p+2j$')
+ title('$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$')
+ subplot(2,1,2)
+ imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ji)));
+ set(gca,'YDir','normal')
+% pclr = pcolor(XX,YY,plt(Ne_ST));
+% set(pclr, 'edgecolor','none'); hold on;
+ xlabel('$t$'); ylabel('$p+2j$')
+ title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$')
+ suptitle(DATA.param_title);
+
+end
+
+end
+
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
new file mode 100644
index 0000000000000000000000000000000000000000..d57f6755ee5dab68e31e0e3689c51d1453a72a64
--- /dev/null
+++ b/matlab/plot/spectrum_1D.m
@@ -0,0 +1,129 @@
+function [ FIGURE ] = spectrum_1D( data, options )
+
+options.PLAN = 'kxky';
+options.COMP = 1;
+options.POLARPLOT = 0;
+toplot = process_field(data,options);
+t = data.Ts3D; frames = toplot.FRAMES;
+
+colors = jet(numel(frames));
+
+switch options.NAME
+ case '\Gamma_x'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS];
+ yname = '$\sum_{k_y}|\Gamma_k|$'; fieldname = 'transport';
+ case '\phi'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS];
+ yname = '$\sum_{k_y}|\phi|$'; fieldname = 'ES pot.';
+ case 'n_i'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS];
+ yname = '$\sum_{k_y}|n_i|$'; fieldname = 'ion dens.';
+ case 'n_e'
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS];
+ yname = '$\sum_{k_y}|n_e|$'; fieldname = 'elec. dens.';
+end
+
+switch options.COMPXY
+ case 'avg'
+ compx = @(x) mean(x,1);
+ compy = @(x) mean(x,2);
+ case 'sum'
+ compx = @(x) sum(x,1);
+ compy = @(x) sum(x,2);
+ case 'max'
+ compx = @(x) max(x,1);
+ compy = @(x) max(x,2);
+ otherwise
+ compx = @(x) x(data.kx==0,:);
+ compy = @(x) x(:,data.ky==0);
+end
+
+set(gcf, 'Position', [20 50 5000 2000])
+subplot(1,2,1)
+
+ k = data.kx;
+ xname = '$k_x$';
+
+ nmax = ceil(data.Nkx*2/3);
+ shiftx = @(x) x;%(1:nmax);
+ shifty = @(x) x;%(1:nmax);
+ switch options.COMPT
+ case 'avg'
+ it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
+ Gk = compy(abs(mean(toplot.FIELD(:,:,it0:it1),3)));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ X = shiftx(k);
+ if options.OK
+ Y = shifty(Gk)./X;
+ else
+ Y = shifty(Gk);
+ end
+ plot(X,Y,'DisplayName','t-averaged')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Gk = compy(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ X = shiftx(k);
+ if options.OK
+ Y = shifty(Gk)./X;
+ else
+ Y = shifty(Gk);
+ end
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ end
+ grid on
+ title(['HeLaZ $k_x$ ',fieldname,' spectrum']); legend('show','Location','eastoutside')
+ xlabel(xname); ylabel(yname)
+
+subplot(1,2,2)
+
+ k = data.ky;
+ xname = '$k_y$';
+ nmax = floor(data.Nky/2*2/3);
+ AA = @(x) x(1:nmax);
+ shiftx = @(x) x;%AA(x);
+ shifty = @(x) x;%AA(ifftshift(x));
+ switch options.COMPT
+ case 'avg'
+ it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
+ Gk = compx(abs(mean(toplot.FIELD(:,:,it0:it1),3)))';
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ X = k(k>0); X = X(1:end-1);
+ Yp= Gk(k>0);
+ Ym= Gk(k<0);
+ Y = Yp(1:end-1) + Ym(end:-1:1);
+ Y = Y(end:-1:1);
+ plot(X,Y,'DisplayName','t-averaged')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Gk = compx(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ X = k(k>0); X = X(1:end-1);
+ Yp= Gk(k>0);
+ Ym= Gk(k<0);
+ Y = Yp(1:end-1) + Ym(end:-1:1);
+ Y = Y(end:-1:1);
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ end
+ grid on
+% title('HeLaZ $k_y$ transport spectrum');
+ legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(yname)
+end
+
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
new file mode 100644
index 0000000000000000000000000000000000000000..5a9d0653fc038f03433cc666b1be230a5b529d74
--- /dev/null
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -0,0 +1,34 @@
+function [ fig ] = statistical_transport_averaging( data, options )
+scale = (1/data.Nx/data.Ny)^2;
+Trange = options.T;
+[~,it0] = min(abs(Trange(1)-data.Ts0D));
+[~,it1] = min(abs(Trange(end)-data.Ts0D));
+gamma = data.PGAMMA_RI(it0:it1)*scale;
+dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
+if ~dt_const
+ disp('DT not const on given interval');
+else
+
+ Ntot = (it1-it0)+1;
+
+ transp_seg_avg = 1:Ntot;
+ transp_seg_std = 1:Ntot;
+
+ for Np = 1:Ntot % Loop on the number of segments
+ transp_seg_avg(Np) = mean(gamma(1:Np));
+ transp_seg_std(Np) = std(gamma(1:Np));
+ end
+
+ time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0));
+
+ fig = figure;
+% subplot(211)
+ plot(time_seg,transp_seg_avg,'-'); hold on;
+ xlabel('Averaging time'); ylabel('transport value');
+
+% subplot(212)
+% errorbar(N_seg,transp_seg_avg,transp_seg_std);
+% xlabel('Averaging #points'); ylabel('transport value');
+end
+end
+
diff --git a/matlab/plot/transport_spectrum.m b/matlab/plot/transport_spectrum.m
new file mode 100644
index 0000000000000000000000000000000000000000..76812d635c7abd5e4982ce85cd3f4744a83f975e
--- /dev/null
+++ b/matlab/plot/transport_spectrum.m
@@ -0,0 +1,63 @@
+function [ FIGURE ] = transport_spectrum( data, options )
+
+options.NAME = '\Gamma_x';
+% options.NAME = '\phi';
+% options.NAME = 'n_i';
+options.PLAN = 'kxky';
+options.COMP = 1;
+
+toplot = process_field(data,options);
+t = data.Ts3D; frames = toplot.FRAMES;
+
+colors = jet(numel(frames));
+
+FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS];
+set(gcf, 'Position', [20 50 5000 2000])
+subplot(1,2,1)
+
+ sdim = 2;
+ k = data.kx;
+ xname = '$k_x$';
+ yname = '$\sum_{k_y}|\Gamma_k|$';
+ nmax = ceil(data.Nkx*2/3);
+ shiftx = @(x) x(1:nmax);
+ shifty = @(x) x(1:nmax);
+
+ for it = 1:numel(toplot.FRAMES)
+ Gk = sum(abs(toplot.FIELD(:,:,toplot.FRAMES(it))),sdim);
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ semilogy(shiftx(k), shifty(Gk),'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ grid on
+ title('HeLaZ $k_x$ transport spectrum'); legend('show','Location','eastoutside')
+ xlabel(xname); ylabel(yname)
+
+subplot(1,2,2)
+
+ sdim = 1;
+ k = data.ky;
+ xname = '$k_y$';
+ yname = '$\sum_{k_x}|\Gamma_k|$';
+ nmax = floor(data.Nky/2*2/3);
+ AA = @(x) x(1:nmax);
+ shiftx = @(x) AA(x);
+ shifty = @(x) AA(ifftshift(x));
+
+ for it = 1:numel(toplot.FRAMES)
+ Gk = sum(abs(toplot.FIELD(:,:,toplot.FRAMES(it))),sdim);
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(abs(Gk));
+ end
+ semilogy(shiftx(k), shifty(Gk),'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ grid on
+ title('HeLaZ $k_y$ transport spectrum'); legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(yname)
+end
+
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 376186a7c062ecea479485046544d2acdb1fd679..ff8d00eeb4d915741bcfef91e9e0c85eee7c163a 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -38,6 +38,10 @@ switch OPTIONS.PLAN
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
+ case 'sx'
+ XNAME = '$v_\parallel$'; YNAME = '$\mu$';
+ [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
+ REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0;
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
@@ -316,34 +320,82 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case 's_y'
+ NAME = 'sy';
+ [~, KX] = meshgrid(DATA.ky, DATA.kx);
+ vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ for it = FRAMES % Compute the 3D real transport for each frame
+ for iz = 1:DATA.Nz
+ vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
+ end
+ end
+ if REALP % plot in real space
+ for it = FRAMES
+ FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,it))));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
case '\Gamma_x'
- NAME = 'Gx';
- [KY, ~] = meshgrid(DATA.ky, DATA.kx);
- Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
- for it = FRAMES % Compute the 3D real transport for each frame
- for iz = 1:DATA.Nz
- Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))...
- .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ NAME = 'Gx';
+ [KY, ~] = meshgrid(DATA.ky, DATA.kx);
+ Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ for it = FRAMES % Compute the 3D real transport for each frame
+ for iz = 1:DATA.Nz
+ Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))...
+ .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
end
- end
- if REALP % plot in real space
+ if REALP % plot in real space
+ for it = FRAMES
+ FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ shift_x = @(x) fftshift(x,2);
+ shift_y = @(x) fftshift(x,2);
+ tmp = zeros(DATA.Nx,DATA.Ny,Nz);
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
+ end
+ end
+ case 'f_i'
+ shift_x = @(x) x; shift_y = @(y) y;
+ NAME = 'fi'; OPTIONS.SPECIE = 'i';
+ [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
+ [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
+ dit_ = 1;%ceil((it1_-it0_+1)/10);
+ FRAMES = it0_:dit_:it1_;
+ iz = 1;
for it = FRAMES
- FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
- tmp = zeros(DATA.Nx,DATA.Ny,Nz);
+ OPTIONS.T = DATA.Ts5D(it);
+ [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
+ end
+ case 'f_e'
+ shift_x = @(x) x; shift_y = @(y) y;
+ NAME = 'fe'; OPTIONS.SPECIE = 'e';
+ [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
+ [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
+ dit_ = 1;%ceil((it1_-it0_+1)/10);
+ FRAMES = it0_:dit_:it1_;
+ iz = 1;
for it = FRAMES
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
+ OPTIONS.T = DATA.Ts5D(it);
+ [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
- end
end
TOPLOT.FIELD = FIELD;
+TOPLOT.FRAMES = FRAMES;
TOPLOT.X = shift_x(X);
TOPLOT.Y = shift_y(Y);
TOPLOT.FIELDNAME = FIELDNAME;
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 5323138fb0a971546dd1c31467e9ce96287dcbb7..6bb65b00fd6d776d7788d9c17bc9fa95771f4c54 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -62,7 +62,7 @@ MODULE grid
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kyarray, kyarray_full
! Kperp array depends on kx, ky, z (geometry), eo (even or odd zgrid)
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC :: kparray
- REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max!, kp_max
+ REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max
REAL(dp), PUBLIC, PROTECTED :: local_kxmax, local_kymax
INTEGER, PUBLIC, PROTECTED :: ikxs, ikxe, ikys, ikye!, ikps, ikpe
INTEGER, PUBLIC, PROTECTED :: ikx_0, iky_0, ikx_max, iky_max ! Indices of k-grid origin and max
@@ -242,9 +242,11 @@ CONTAINS
IF (Nx .EQ. 1) THEN
deltakx = 0._dp
kx_max = 0._dp
+ kx_min = 0._dp
ELSE
deltakx = 2._dp*PI/Lx
kx_max = Nkx*deltakx
+ kx_min = deltakx
ENDIF
! Build the full grids on process 0 to diagnose it without comm
@@ -313,11 +315,13 @@ CONTAINS
contains_ky0 = .true.
ky_max = 0._dp
iky_max = 1
+ ky_min = 0._dp
kyarray_full(1) = 0._dp
local_kymax = 0._dp
ELSE ! Build apprpopriate grid
deltaky = 2._dp*PI/Ly
ky_max = (Ny/2)*deltakx
+ ky_min = deltaky
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
local_kymax = 0._dp
DO iky = ikys,ikye
diff --git a/src/inital.F90 b/src/inital.F90
index 06a336b0f6ebecc7e11c6727d27a837ea70e7d6a..c103b9149fe3a8c5eb01dec4e02db36082924224 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -46,6 +46,10 @@ SUBROUTINE inital
IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
CALL initialize_blob
CALL poisson ! get phi_0 = phi(N_0)
+ CASE('ppj')
+ IF (my_id .EQ. 0) WRITE(*,*) 'ppj init ~ GENE'
+ call init_ppj
+ CALL poisson ! get phi_0 = phi(N_0)
END SELECT
ENDIF
@@ -356,3 +360,124 @@ SUBROUTINE initialize_blob
ENDDO
END SUBROUTINE initialize_blob
!******************************************************************************!
+
+
+!******************************************************************************!
+!!!!!!! Initialize the gyrocenter in a ppj gene manner (power law)
+!******************************************************************************!
+SUBROUTINE init_ppj
+ USE basic
+ USE grid
+ USE fields
+ USE prec_const
+ USE utility, ONLY: checkfield
+ USE initial_par
+ USE model, ONLY : LINEARITY, KIN_E
+ IMPLICIT NONE
+
+ REAL(dp) :: noise
+ REAL(dp) :: kx, ky, sigma, gain, ky_shift
+ INTEGER, DIMENSION(12) :: iseedarr
+
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
+
+ !**** Broad noise initialization *******************************************
+ ! Electrons
+ IF (KIN_E) THEN
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+ DO ikx=ikxs,ikxe
+ kx = kxarray(ikx)
+ DO iky=ikys,ikye
+ ky = kyarray(iky)
+ DO iz=izs,ize
+ IF (kx .NE. 0) THEN
+ IF(ky .NE. 0) THEN
+ moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
+ ELSE
+ moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
+ ENDIF
+ ELSE
+ IF(ky .NE. 0) THEN
+ moments_e(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ ELSE
+ moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
+ ENDIF
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ IF ( contains_kx0 ) THEN
+ DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+ moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+ END DO
+ ENDIF
+ ELSE
+ moments_e(ip,ij,:,:,:,:) = 0._dp
+ ENDIF
+ END DO
+ END DO
+ ENDIF
+
+ ! Ions
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+ DO ikx=ikxs,ikxe
+ kx = kxarray(ikx)
+ DO iky=ikys,ikye
+ ky = kyarray(iky)
+ DO iz=izs,ize
+ IF (kx .NE. 0) THEN
+ IF(ky .NE. 0) THEN
+ moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
+ ELSE
+ moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
+ ENDIF
+ ELSE
+ IF(ky .NE. 0) THEN
+ moments_i(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ ELSE
+ moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
+ ENDIF
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ IF ( contains_kx0 ) THEN
+ DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+ moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+ END DO
+ ENDIF
+ ELSE
+ moments_i(ip,ij,:,:,:,:) = 0._dp
+ ENDIF
+ END DO
+ END DO
+
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz=izs,ize
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+END SUBROUTINE init_ppj
+!******************************************************************************!
diff --git a/wk/Dimits_shift_results.m b/wk/Dimits_shift_results.m
index 5bbd02a31c811713a48a5e41f9076be9b0de41a3..3f846c777f1175e9dcb78185f4281f759fadeb8a 100644
--- a/wk/Dimits_shift_results.m
+++ b/wk/Dimits_shift_results.m
@@ -1,24 +1,24 @@
%% Nu = 1e-2
-KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LD = [5.5e-2 1.3e-1 2.7e-1 4.6e-1 1.0e+0];
-conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.5];
+Ginf_LDGK = [5.5e-2 1.3e-1 2.7e-1 4.6e-1 1.0e+0 0.0e+0 0.0e+0];
+conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_SG = [1.2e-3 2.8e-3 1.7e-2 8.9e-2 3.0e+0];
-conv_SG = [9.0e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_SGGK = [1.2e-3 2.8e-3 1.7e-2 8.9e-2 3.0e+0 0.0e+0 0.0e+0];
+conv_SGGK = [9.0e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_DG = [1.1e-4 5.6e-4 9.0e-4 6.5e-3 7.0e+0];
-conv_DG = [1.6e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_DGGK = [1.1e-4 5.6e-4 9.0e-4 6.5e-3 7.0e+0 0.0e+0 0.0e+0];
+conv_DG = [1.6e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
-semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,conv_SG,'x-k','DisplayName','conv'); hold on;
+semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,conv_SGGK,'x-k','DisplayName','conv'); hold on;
-semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
semilogy(KN,conv_DG,'x-k','DisplayName','conv'); hold on;
title('$\nu=0.01$');
@@ -27,34 +27,41 @@ ylabel('Transport');
%% Nu = 5e-2
KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_SG = [0.0e+0 0.0e+0 0.0e+0 2.0e-1 0.0e+0];
-conv_SG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_DG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_LDGK = [0.0e+0 0.0e+0 0.0e+0 8.3e-1 0.0e+0];
+conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+Ginf_SGGK = [0.0e+0 0.0e+0 0.0e+0 2.0e-1 0.0e+0];
+conv_SGGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+Ginf_SGDK = [0.0e+0 0.0e+0 0.0e+0 5.6e-3 0.0e+0];
+conv_SGDK = [0.0e+0 0.0e+0 0.0e+0 2.0e-3 0.0e+0];
+
+Ginf_DGGK = [0.0e+0 0.0e+0 0.0e+0 1.5e-2 0.0e+0];
figure
-semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
-semilogy(KN,conv_SG,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau GK'); hold on;
+semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama GK'); hold on;
+semilogy(KN,Ginf_SGDK,'o-c','DisplayName','Sugama DK'); hold on;
+semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty GK'); hold on;
+semilogy(KN,conv_LDGK,'x--g','DisplayName','changing params'); hold on;
+semilogy(KN,conv_SGGK,'x--b','DisplayName','changing params'); hold on;
+semilogy(KN,conv_SGDK,'x--c','DisplayName','changing params'); hold on;
title('$\nu=0.05$');
xlabel('Drive');
ylabel('Transport');
%% Nu = 1e-1
-KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LD = [2.5e-2 2.5e-1 6.3e-1 1.3e+0 2.3e+0];
-conv_LD = [0.0e+0 2.5e-1 0.0e+0 1.2e+0 0.0e+0];
-Ginf_SG = [0.0e-9 1.3e-2 9.1e-2 3.0e-1 7.8e-1];
-Ginf_DG = [2.0e-4 6.4e-3 3.6e-2 3.1e-1 0.0e+0];
+KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.5];
+Ginf_LDGK = [2.5e-2 2.5e-1 6.3e-1 1.3e+0 2.3e+0 4.3e+0 0.0e+0];
+conv_LDGK = [0.0e+0 2.5e-1 0.0e+0 1.2e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_SGGK = [0.0e-9 1.3e-2 9.1e-2 3.0e-1 7.8e-1 1.4e+0 3.5e+1];
+Ginf_DGGK = [2.0e-4 6.4e-3 3.6e-2 3.1e-1 0.0e+0 3.0e-1 3.3e+1];
figure
-semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_LDGK,'x--k','DisplayName','changing params'); hold on;
title('$\nu=0.1$');
xlabel('Drive');
ylabel('Transport');
@@ -62,18 +69,18 @@ ylabel('Transport');
%% Nu = 5e-1
KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_SG = [0.0e-9 3.9e-2 2.7e-1 6.6e-1 1.6e+0];
-conv_SG = [0.0e-9 0.0e+0 3.0e-1 0.0e+0 1.6e+0];
-Ginf_DG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_SGGK = [0.0e-9 3.9e-2 2.7e-1 6.6e-1 1.6e+0];
+conv_SGGK = [0.0e-9 0.0e+0 3.0e-1 0.0e+0 1.6e+0];
+Ginf_DGGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
-semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
-semilogy(KN,conv_SG,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_LDGK,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,conv_SGGK,'x--k','DisplayName','changing params'); hold on;
title('$\nu=0.5$');
xlabel('Drive');
ylabel('Transport');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 941a4ace3d2972ed825b7350413e982efe8e47c5..d3e9d353f6a6a3698b0e1efce5c24683851163a4 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,109 +1,14 @@
addpath(genpath('../matlab')) % ... add
-addpath(genpath('../matlab/plots')) % ... add
-outfile ='';
-%% Directory of the simulation
-if 1% Local results
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-%% nu = 5e-1
-% Sugama
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK';
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK';
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4
-
-%% nu = 1e-1
-% Landau
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK';
-
-% Sugama
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK';
-
-% Dougherty
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK';
-
-%% nu = 5e-2
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD)
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
-
-% testing various NL closures
-% outfile ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
-outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-
-%% nu = 1e-2
-% Landau
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK';
-
-% Sugama
-% outfile ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD)
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
-% Dougherty
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK';
-
-%% nu = 5e-3
-% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
-% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
-% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
-% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
-
-%% nu = 0
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120x240_kN_1.6_kT_0.4_nu_0_mu_0.05';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_0_mu_0.05';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_5e-2';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_11x6_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
-
- LOCALDIR = ['../results/',outfile,'/'];
- MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
- system(['mkdir -p ',MISCDIR]);
- CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
- system(CMD);
-else% Marconi results
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
-% BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
-MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
-end
+addpath(genpath('../matlab/plot')) % ... add
+addpath(genpath('../matlab/compute')) % ... add
+addpath(genpath('../matlab/load')) % ... add
%% Load the results
+LOCALDIR = ['../results/',outfile,'/'];
+MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
+system(['mkdir -p ',MISCDIR]);
+CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
+system(CMD);
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -116,18 +21,16 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 =100; TAVG_1 = 80000; % Averaging times duration
+TAVG_0 = 500; TAVG_1 = 10000; % Averaging times duration
% chose your field to plot in spacetime diag (uzf,szf,Gx)
-fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi');
+fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'v_y',1);
save_figure(data,fig)
end
if 0
-%% ZF caracteristics (space time diagrams
-TAVG_0 = 500; TAVG_1 = 18000; % Averaging times duration
-% chose your field to plot in spacetime diag (uzf,szf,Gx)
-fig = ZF_spacetime(data,TAVG_0,TAVG_1);
-save_figure(data,fig)
+%% statistical transport averaging
+options.T = [1500 5000];
+fig = statistical_transport_averaging(data,options);
end
if 0
@@ -141,9 +44,11 @@ options.NAME = '\phi';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xy';
+% options.NAME = 'f_e';
+% options.PLAN = 'sx';
options.COMP = 1;
-options.TIME =200:3:800;
-% options.TIME = 140:0.5:160;
+% options.TIME = data.Ts5D;
+options.TIME = 0:3:1000;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
end
@@ -153,14 +58,17 @@ if 0
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
-options.AXISEQUAL = 0;
+options.AXISEQUAL = 1;
% options.NAME = '\phi';
+% options.NAME = 'v_y';
+options.NAME = 'n_e^{NZ}';
% options.NAME = '\Gamma_x';
-options.NAME = 'n_i^{NZ}';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
+% options.NAME = 'f_e';
+% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = 1.2e4+[600 800 1000];
+options.TIME = [550:20:800];
data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -171,136 +79,85 @@ if 0
options.INTERP = 1;
options.NAME = 'n_i';
options.PLANES = 1;
-options.TIME = 1.2e4+[0 500 1000];
+options.TIME = [0 500 1000];
data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
save_figure(data,fig)
end
if 0
-%%
-TAVG_0 = 1000; TAVG_1 = 5000; % Averaging times duration
-ZF_fourier_analysis
+%% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
+options.SPAR = linspace(-3,3,64)+(6/127/2);
+options.XPERP = linspace( 0,6,64);
+% options.SPAR = vp';
+% options.XPERP = mu';
+options.Z = 1;
+options.T = 4000;
+options.CTR = 1;
+options.ONED = 0;
+fig = plot_fa(data,options);
+save_figure(data,fig)
end
if 0
-%% Radial shear profile (with moving average)
-tf = 1000+[0:100:1000];
-ymax = 0;
-figure
-for i_ = 1:numel(tf)
-[~,it] = min(abs(Ts3D-tf(i_)));
-data = squeeze((mean(dx2phi(:,:,1,it),2)));
-step = 50;
-plot(movmean(x,step),movmean(data,step),'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on;
-ymax = max([ymax abs(min(data)) abs(max(data))]);
-end
-xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax);
-xlabel('$x/\rho_s$'); ylabel('$s_{E\times B,x}$'); grid on
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.TIME = 1000:4000;
+options.P2J = 1;
+options.ST = 1;
+options.NORMALIZED = 0;
+fig = show_moments_spectrum(data,options);
+save_figure(data,fig)
end
if 0
-%% zonal vs nonzonal energies for phi(t)
-t0 = 200; [~, it0] = min(abs(Ts3D-t0));
-itend = Ns3D;
-trange = it0:itend;
-pltx = @(x) x;%-x(1);
-plty = @(x) x./max(squeeze(x));
-fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS];
-set(gcf, 'Position', [100, 100, 1400, 500])
-subplot(121)
-% yyaxis left
- semilogy(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on;
-% yyaxis right
- semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName',['NZ, $k_p>0$, ',CONAME]);
- semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName',['NZ, $k_p>1$, ',CONAME]);
- semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName',['NZ, $k_p>2$, ',CONAME]);
-% semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]);
- title('Energy'); legend('Location','southeast')
- xlim([Ts3D(it0), Ts3D(itend)]);
- ylim([1e-3, 1.5])
- xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
-
-subplot(122)
-% plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)));
- plot3(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)),Ts3D(trange));
- title('Phase space'); legend(CONAME)
- xlabel('$E_Z$'); ylabel('$E_{NZ}$'); zlabel('time'); grid on;% xlim([0 500]);
+%% 1D spectrum
+options.TIME = 1000:3000;
+options.NORM = 1;
+% options.NAME = '\phi';
+% options.NAME = 'n_i';
+options.NAME ='\Gamma_x';
+options.PLAN = 'kxky';
+options.COMPZ = 1;
+options.OK = 0;
+options.COMPXY = 'max';
+options.COMPT = 'avg';
+options.PLOT = 'semilogy';
+fig = spectrum_1D(data,options);
+% save_figure(data,fig)
end
if 0
-%% Conservation laws
-Nxmax = Nx/2;
-Nymax = Ny/2;
-mflux_x_i = squeeze(sum((Gamma_x( 1,1:Nxmax,:)+Gamma_x( 1,2:Nxmax+1,:))/2,2)./sum(Gamma_x( 1,1:Nxmax,:)));
-mflux_x_o = squeeze(sum((Gamma_x( Nxmax,1:Nxmax,:)+Gamma_x( Nxmax,2:Nxmax+1,:))/2,2)./sum(Gamma_x( Nxmax,1:Nxmax,:)));
-mflux_y_i = squeeze(sum((Gamma_y(1:Nxmax, 1,:)+Gamma_y(2:Nxmax+1, 1,:))/2,1)./sum(Gamma_y(1:Nxmax, 1,:)));
-mflux_y_o = squeeze(sum((Gamma_y(1:Nxmax, Nymax,:)+Gamma_y(2:Nxmax+1, Nymax,:))/2,1)./sum(Gamma_y(1:Nxmax, Nymax,:)));
-
-mass_cons = mflux_x_i - mflux_x_o + mflux_y_i - mflux_y_o;
-%%
-figure
-plt = @(x) squeeze(mean(mean(x(:,:,1,:),1),2));
-subplot(211)
- plot(Ts3D,plt(dens_e+dens_i),'DisplayName','$\delta n_e + \delta n_i$'); hold on;
- plot(Ts3D,plt(ne00+ni00),'DisplayName','$\delta n_e^{00} + \delta n_i^{00}$'); hold on;
- plot(Ts3D,plt(temp_e+temp_i),'DisplayName','$\delta T_e + \delta T_i$'); hold on;
- legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)])
-subplot(212);
- plot(Ts3D,mass_cons*(2*pi/Nx/Ny)^2,'DisplayName','in-out'); hold on
-% plot(Ts3D,squeeze(mflux_x_i),'DisplayName','Flux i x');
-% plot(Ts3D,squeeze(mflux_x_o),'DisplayName','Flux o x');
-% plot(Ts3D,squeeze(mflux_y_i),'DisplayName','Flux i y');
-% plot(Ts3D,squeeze(mflux_y_o),'DisplayName','Flux o y');
- legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]); %ylim([-0.1, 2]*mean(mflux_x_i))
+%% 1D real plot
+options.TIME = 1500:100:2500;
+options.NORM = 0;
+% options.NAME = '\phi';
+% options.NAME = 'n_i';
+% options.NAME ='\Gamma_x';
+options.NAME ='s_y';
+options.COMPZ = 1;
+options.COMPXY = 'avg';
+options.COMPT = 'avg';
+options.MOVMT = 1;
+fig = real_plot_1D(data,options);
+% save_figure(data,fig)
end
if 0
-%% Zonal profiles (ky=0)
-
-% Chose the field to plot
-FIELD = data.Ne00.*conj(data.Ne00); FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2'
-% FIELD = Ni00.*conj(Ni00); FNAME = 'Ni002'; FIELDLTX = '|N_i^{00}|^2'
-% FIELD = abs(PHI); FNAME = 'absPHI'; FIELDLTX = '|\tilde\phi|';
-% FIELD = PHI.*conj(PHI); FNAME = 'PHI2'; FIELDLTX = '|\tilde\phi|^2';
-% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
-
-% Chose when to plot it
-tf = 200:200:1200;
-% tf = 8000;
-
-% Sliced
-plt = @(x,it) x( 2:end, 1,1,it)./max(max(x( 2:end, 1,1,:))); XNAME = 'k_x';
-%
-TNAME = [];
-fig = figure; FIGNAME = ['Zonal_',FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 600,400])
-plt_2 = @(x) (fftshift(x,2));
- for i_ = 1:numel(tf)
- [~,it] = min(abs(data.Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(data.Ts3D(it))];
- d2plot = plt_2(squeeze(plt(FIELD,it)));
- semilogy(data.kx(2:end),d2plot,'-','DisplayName',sprintf('$t c_s/R=%.0f$',data.Ts3D(it))); hold on; grid on;
- xlabel(latexize(XNAME));
- end
-title(['$',FIELDLTX,'$ Zonal Spectrum']); legend('show');
-
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 0.01:0.01:1.0;
+options.TIME =500:1:3000;
+options.NMA = 1;
+options.NMODES = 20;
+fig = mode_growth_meter(data,options);
+save_figure(data,fig)
end
if 0
-%% Time evolutions and growth rate
-plot_time_evolution_and_gr
+%% ZF caracteristics (space time diagrams
+TAVG_0 = 200; TAVG_1 = 3000; % Averaging times duration
+% chose your field to plot in spacetime diag (uzf,szf,Gx)
+fig = ZF_spacetime(data,TAVG_0,TAVG_1);
+save_figure(data,fig)
end
-
-if 0
-%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 0000; tend = 1000;
-% tstart = 10000; tend = 12000;
-% Chose the field to plot
-% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
-% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0;
-plot_kperp_spectrum
-end
\ No newline at end of file
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
new file mode 100644
index 0000000000000000000000000000000000000000..aa87383b5a00f19bb9efbaa6606c92856d59b178
--- /dev/null
+++ b/wk/analysis_header.m
@@ -0,0 +1,154 @@
+%% Directory of the simulation
+% if 1% Local results
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+% outfile ='debug/ppj_init';
+%% nu = 5e-1
+% Sugama
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_7x4_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4
+
+%% nu = 1e-1
+% Landau
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK';
+
+% Sugama
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK';
+
+% Dougherty
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK';
+
+%% nu = 5e-2
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD)
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+
+% testing various NL closures
+% outfile ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+
+% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+
+% outfile ='Hallenbert_nu_5e-02/256x64_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/256x64_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_21x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_17x9_L_120_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+
+% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
+% outfile ='Hallenbert_nu_5e-02/128x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_FCGK';
+
+%% nu = 1e-2
+% Landau
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK';
+
+% Sugama
+% outfile ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD)
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% Dougherty
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK';
+
+%% nu = 5e-3
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+
+%% nu = 0
+
+outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLDGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLRGK_0.1_muxy_1e-2';
+
+
+% outfile ='Hallenbert_fig2b/200x32_11x6_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuSGGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuDGGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLDGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLRGK_0.1_muxy_1e-1';
+
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLDGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_9x5_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
+
+%% Transport scan
+% outfile = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5';
+
+
+% outfile = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1';
+% outfile = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/SG_10x5_conv_test';
+% outfile = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_Lx1.5';
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_HD';
+
+% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
+
+% outfile = 'predator_prey/SG_Kn_1.7_nu_0.01';
+
+% else% Marconi results
+% outfile ='';
+% outfile ='';
+% outfile ='';
+% outfile ='';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt';
+% % outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+% % BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
+% MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
+% end
+
+analysis_3D
\ No newline at end of file
diff --git a/wk/benchmark_HeLaZ_gene_transport_zpinch.m b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
new file mode 100644
index 0000000000000000000000000000000000000000..0359ffb73741d5073bf2fc86533252e6d25714e3
--- /dev/null
+++ b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
@@ -0,0 +1,131 @@
+%% Nu = 0
+KN = [ 1.6 2.0 2.5];
+Ginf_200x32x05x03 = [2.8e-3 3.1e-1 2.4e+1];
+Gstd_200x32x05x03 = [4.4e-3 2.9e-1 2.0e+0];
+Gmin_200x32x05x03 = [6.1e-5 2.5e-3 1.9e+1];
+Gmax_200x32x05x03 = [2.3e-2 1.2e+0 2.9e+1];
+
+Ginf_200x32x11x06 = [7.3e-3 6.3e-1 2.4e+1];
+Gstd_200x32x11x06 = [1.3e-2 4.0e-1 5.9e+0];
+Gmin_200x32x11x06 = [1.8e-5 4.8e-3 1.7e+1];
+Gmax_200x32x11x06 = [1.7e-2 2.0e+0 3.1e+1];
+
+Ginf_200x32x21x11 = [1.2e-2 0 0];
+
+Ginf_GENE = [8.5e-3 5.5e-1 2.3e+1];
+Gstd_GENE = [5.1e-3 3.7e-1 2.6e+0];
+Gmin_GENE = [8.0e-4 9.2e-3 2.0e+1];
+Gmax_GENE = [6.0e-3 1.5e+0 2.3e+1];
+
+
+fig = figure; set(gcf,'Position',[250 250 600 300]);
+% plot(KN,Ginf_200x32x05x03./Ginf_GENE,'v','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% plot(KN,Ginf_200x32x11x06./Ginf_GENE,'^','DisplayName','HeLaZ 200x32x10x05','MarkerSize',10); hold on;
+% errorbar(KN-0.05,Ginf_200x32x05x03,Gstd_200x32x11x06/2,'v','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% errorbar(KN+0.05,Ginf_200x32x11x06,Gstd_200x32x11x06/2,'v','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% errorbar(KN ,Ginf_GENE ,Gstd_GENE/2 ,'xk','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% errorbar(KN-0.05,Ginf_200x32x05x03,Gmin_200x32x05x03,Gmax_200x32x05x03,'v','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% errorbar(KN+0.05,Ginf_200x32x11x06,Gmin_200x32x11x06,Gmax_200x32x11x06,'v','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% errorbar(KN ,Ginf_GENE ,Gmin_GENE ,Gmax_GENE ,'xk','DisplayName','HeLaZ 200x32x05x03','MarkerSize',10); hold on;
+% xlabel('$\kappa_N$');
+% ylabel('$\Gamma_x^\infty$');
+
+loglog(Ginf_GENE,Ginf_200x32x05x03,'*','DisplayName','HeLaZ 5x3','MarkerSize',10); hold on;
+loglog(Ginf_GENE,Ginf_200x32x11x06,'*','DisplayName','HeLaZ 10x5','MarkerSize',8); hold on;
+loglog(Ginf_GENE,Ginf_200x32x21x11,'*','DisplayName','HeLaZ 20x10','MarkerSize',8); hold on;
+% title('Collisionless transport benchmark');
+xy = [1e-3 1e+2];
+loglog(xy,xy,'--k','DisplayName','x=y');
+% loglog(Ginf_GENE,Gmin_GENE,'-k','DisplayName','x=y');
+% loglog(Ginf_GENE,Gmax_GENE,'-k','DisplayName','x=y');
+xlabel('Gene $\Gamma_x^\infty$');
+ylabel('HeLaZ $\Gamma_x^\infty$');
+xlim(xy); ylim(xy);
+grid on
+% saveas(fig,'/home/ahoffman/Dropbox/Applications/Overleaf/Paper 1/figures/transport_benchmark.eps')
+
+
+%% Nu = 0.1
+fig = figure; set(gcf,'Position',[250 250 600 300]);
+
+KN = [ 1.6 2.0 2.5];
+LDGK_200x32x05x03 = [2.6e-1 4.3e+0 3.8e+1];
+LRGK_200x32x05x03 = [4.0e-1 2.4e+0 2.9e+1];
+SGGK_200x32x05x03 = [1.4e-2 1.4e+0 3.5e+1];
+DGGK_200x32x05x03 = [3.8e-3 1.8e-1 3.0e+1];
+
+SGGK_200x32x11x06 = [1.5e-2 1.3e+0 0.0e+0];
+DGGK_200x32x11x06 = [2.4e-3 3.8e-1 0.0e+0];
+
+X = Ginf_GENE;
+% X = KN;
+clr_ = line_colors;
+
+loglog(X,DGGK_200x32x05x03,'v','DisplayName','Dougherty','MarkerSize',10,'Color',clr_(2,:)); hold on;
+loglog(X,SGGK_200x32x05x03,'s','DisplayName', 'Sugama','MarkerSize',10,'Color',clr_(1,:)); hold on;
+loglog(X,LDGK_200x32x05x03,'d','DisplayName', 'Landau','MarkerSize',10,'Color',clr_(5,:)); hold on;
+loglog(X,LRGK_200x32x05x03,'^','DisplayName', 'Lorentz','MarkerSize',10,'Color',clr_(3,:)); hold on;
+
+loglog(X,DGGK_200x32x11x06,'vk','DisplayName', 'conv','MarkerSize',8); hold on;
+loglog(X,SGGK_200x32x11x06,'sk','DisplayName', 'conv','MarkerSize',8); hold on;
+
+%% Nu = 0.1
+fig = figure; set(gcf,'Position',[250 250 600 300]);
+
+% KN = [ 1.6 2.0 2.5];
+% LDGK_200x32x05x03 = [2.6e-1 3.6e+0 3.8e+1];
+% LRGK_200x32x05x03 = [4.0e-1 2.0e+0 2.9e+1];
+% SGGK_200x32x05x03 = [1.5e-2 1.4e+0 3.5e+1];
+% DGGK_200x32x05x03 = [3.6e-3 1.8e-1 3.2e+1];
+% NOCO_200x32x05x03 = [3.2e-3 3.1e-1 2.4e+1];
+
+% LRGK_200x32x09x05 = [0.0e+0 2.0e+0 0.0e+0];
+% SGGK_200x32x11x06 = [1.5e-2 1.3e+0 0.0e+0];
+% DGGK_200x32x11x06 = [2.4e-3 3.8e-1 0.0e+0];
+
+% upward scan
+KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
+LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.3e-1 1.3e+0 2.3e+0 3.6e+0 6.3e+0 9.4e+0 1.5e+1 2.3e+1 3.8e+1];
+LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.0e+0 3.5e+0 6.0e+0 4.8e+0 6.0e+0 2.9e+1];
+SGGK_200x32x05x03 = [0.0e-9 1.5e-2 1.0e-1 3.8e-1 8.4e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
+DGGK_200x32x05x03 = [2.0e-4 3.6e-3 1.3e-2 4.0e-2 1.2e-1 1.8e-1 2.0e-1 4.4e-1 6.6e-1 1.3e+1 3.2e+1];
+NOCO_200x32x05x03 = [0.0e+0 1.2e-2 2.6e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 0.0e+0 0.0e+0 0.0e+1 2.4e+1];
+% downward scan
+% KN = [KN 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
+% LDGK_200x32x05x03 = [LDGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+% LRGK_200x32x05x03 = [LRGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+% SGGK_200x32x05x03 = [SGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+% DGGK_200x32x05x03 = [DGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+% NOCO_200x32x05x03 = [NOCO_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 3.9e+0 6.9e+0 1.2e+1 0.0e+0];
+
+
+% X = Ginf_GENE;
+X = KN;
+clr_ = line_colors;
+plt = @(x) x./X;
+% plt = @(x) x;
+msize = 10;
+semilogy(X,plt(DGGK_200x32x05x03),'v','DisplayName','Dougherty','MarkerSize',msize,'Color',clr_(2,:)); hold on;
+semilogy(X,plt(SGGK_200x32x05x03),'s','DisplayName', 'Sugama','MarkerSize',msize,'Color',clr_(1,:)); hold on;
+semilogy(X,plt(LDGK_200x32x05x03),'d','DisplayName', 'Landau','MarkerSize',msize,'Color',clr_(5,:)); hold on;
+semilogy(X,plt(LRGK_200x32x05x03),'^','DisplayName', 'Lorentz','MarkerSize',msize,'Color',clr_(3,:)); hold on;
+semilogy(X,plt(NOCO_200x32x05x03),'*k','DisplayName', '$\nu = 0$','MarkerSize',msize); hold on;
+
+X = [ 1.6 2.0 2.5];
+% loglog(X,DGGK_200x32x11x06./X,'vk','DisplayName', 'conv','MarkerSize',8); hold on;
+% loglog(X,SGGK_200x32x11x06./X,'sk','DisplayName', 'conv','MarkerSize',8); hold on;
+% loglog(X,LRGK_200x32x09x05./X,'^k','DisplayName', 'conv','MarkerSize',8); hold on;
+ylabel('$\Gamma_x^\infty/\kappa_N$, $\nu=0.1$');
+xlabel('$\kappa_N$'); grid on
+xlim([ 1.55 2.55]);
+%%
+% saveas(fig,'/home/ahoffman/Dropbox/Applications/Overleaf/Paper 1/figures/coll_transport_benchmark.eps')
+
+%% Burst study Kn = 1.7
+NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 1.0e-1];
+SGGK_200x32x05x03 = [0.0e+0 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 1.0e-1];
+LRGK_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e-2 0.0e-2 0.0e+0 0.0e+0];
+% NOCO_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+figure
+semilogy(NU,SGGK_200x32x05x03); hold on
+grid on; xlabel('$\nu$'); ylabel('$\Gamma^\infty_x$');
diff --git a/wk/evolve_tracers.m b/wk/evolve_tracers.m
index 05e3613d6657849f92b3987d30bad59c39824a98..2e03ee8e80e1c99a35b7e571cf5c74564151605c 100644
--- a/wk/evolve_tracers.m
+++ b/wk/evolve_tracers.m
@@ -1,12 +1,12 @@
% Options
SHOW_FILM = 0;
-U_TIME =11000; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME =2400; % >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.1;
+dt_ = 0.2;
Nstep = ceil(Tfin/dt_);
% Init tracers
-Np = 40; %number of tracers
+Np = 100; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
@@ -21,12 +21,30 @@ Disp_y = zeros(Np,Nstep);
xmax = max(data.x); xmin = min(data.x);
ymax = max(data.y); ymin = min(data.y);
-% Xp = 0.8*xmax*(0.5-rand(Np));
-% Yp = 0.8*ymax*(0.5-rand(Np));
-dp_ = (xmax-xmin)/(Np-1);
-Xp = linspace(xmin+dp_/2,xmax-dp_/2,Np);
-Yp = zeros(1,Np);
-Zp = zeros(Np);
+INIT = 'round';
+switch INIT
+ case 'lin'
+ % Evenly distributed initial positions
+ dp_ = (xmax-xmin)/(Np-1);
+ Xp = linspace(xmin+dp_/2,xmax-dp_/2,Np);
+ Yp = zeros(1,Np);
+ Zp = zeros(Np);
+ case 'round'
+ % All particles arround a same point
+ xc = 0; yc = 0;
+ theta = rand(1,Np)*2*pi; r = 0.1;
+ Xp = xc + r*cos(theta);
+ Yp = yc + r*sin(theta);
+ Zp = zeros(Np);
+ case 'gauss'
+ % normal distribution arround a point
+ xc = 0; yc = 0; sgm = 1.0;
+ dx = normrnd(0,sgm,[1,Np]); dx = dx - mean(dx);
+ Xp = xc + dx;
+ dy = normrnd(0,sgm,[1,Np]); dy = dy - mean(dy);
+ Yp = yc + dy;
+ Zp = zeros(Np);
+end
% position grid and velocity field
[YY_, XX_ ,ZZ_] = meshgrid(data.y,data.x,data.z);
@@ -41,8 +59,22 @@ ni = zeros(size(XX_));
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));
end
+%% FILM options
+FPS = 30; DELAY = 1/FPS;
+FORMAT = '.gif';
+if SHOW_FILM
+ FILENAME = [data.localdir,'tracer_evolution',FORMAT];
+ switch FORMAT
+ case '.avi'
+ vidfile = VideoWriter(FILENAME,'Uncompressed AVI');
+ vidfile.FrameRate = FPS;
+ open(vidfile);
+ end
+ fig = figure;
+end
%
%Time loop
@@ -50,10 +82,7 @@ t_ = 0;
it = 1;
itu_old = 0;
nbytes = fprintf(2,'frame %d/%d',it,Nstep);
-if SHOW_FILM
- fig = figure;
-end
-while(t_<Tfin)
+while(t_<Tfin && it <= Nstep)
if Evolve_U
[~,itu_] = min(abs(U_TIME+t_-data.Ts3D));
end
@@ -106,7 +135,7 @@ while(t_<Tfin)
iz1 = izC;
end
x0 = data.x(ix0); x1 = data.x(ix1); %left right
- y0 = data.x(iy0); y1 = data.y(iy1); %down top
+ y0 = data.y(iy0); y1 = data.y(iy1); %down top
z0 = data.z(iz0); z1 = data.z(iz1); %back front
if(e_x > 0)
ai__ = (x_ - x0)/(x1-x0); % interp coeff x
@@ -146,8 +175,9 @@ while(t_<Tfin)
u___ = linterp(u__0,u__1,a__i);
- % push the particle
- q = sign(-u___(3));
+% push the particle
+% q = sign(-u___(3));
+ q = -u___(3);
% q =1;
x_ = x_ + dt_*u___(1)*q;
y_ = y_ + dt_*u___(2)*q;
@@ -174,14 +204,14 @@ while(t_<Tfin)
Xp(ip) = x_; Yp(ip) = y_;
end
%% Movie
- if Evolve_U && (itu_old ~= itu_) && SHOW_FILM
+ 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));
scale = max(max(abs(F2P))); % Scaling to normalize
pclr = pcolor(XX_,YY_,F2P/scale);
colormap(bluewhitered);
- set(pclr, 'edgecolor','none'); hold on; caxis([-5,5]);
+ set(pclr, 'edgecolor','none'); hold on; caxis([-2,2]); shading interp
for ip = 1:Np
ia0 = max(1,it-Na);
plot(Traj_x(ip,ia0:it),Traj_y(ip,ia0:it),'.','Color',color(ip,:)); hold on
@@ -194,6 +224,24 @@ while(t_<Tfin)
axis equal
xlim([xmin xmax]); ylim([ymin ymax]);
drawnow
+ % Capture the plot as an image
+ frame = getframe(fig);
+ switch FORMAT
+ case '.gif'
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,32);
+ % Write to the GIF File
+ if it == 1
+ imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ case '.avi'
+ writeVideo(vidfile,frame);
+ otherwise
+ disp('Unknown format');
+ break
+ end
end
t_ = t_ + dt_; it = it + 1; itu_old = itu_;
% terminal info
@@ -203,6 +251,15 @@ while(t_<Tfin)
end
nbytes = fprintf(2,'frame %d/%d',it,Nstep);
end
+disp(' ')
+switch FORMAT
+ case '.gif'
+ disp(['Gif saved @ : ',FILENAME])
+ case '.avi'
+ disp(['Video saved @ : ',FILENAME])
+ close(vidfile);
+end
+
Nt = it;
%% Plot trajectories and statistics
xtot = Disp_x;
@@ -230,11 +287,24 @@ ylabel('$x_p$');
xlim(U_TIME + [0 Tfin]);
subplot(222);
- plot(time_,mean(xtot,1)); hold on
- fit = polyfit(time_,mean(xtot,1),1);
- plot(time_,fit(1)*time_+fit(2),'--k','DisplayName',['$\alpha=',num2str(fit(1)),'$']); hold on
-ylabel('$\langle x \rangle_p$');
-xlim(U_TIME + [0 Tfin]);
+ itf = floor(Nt/2); %fit end time
+ % x^2 displacement
+ plot(time_,mean(xtot.^2,1),'DisplayName','$\langle x.^2\rangle_p$'); hold on
+ fit = polyfit(time_(1:itf),mean(xtot(:,1:itf).^2,1),1);
+ plot(time_,fit(1)*time_+fit(2),'--k'); hold on
+ ylabel('$\langle x^2 \rangle_p$');
+
+% % y^2 displacement
+% fit = polyfit(time_(1:itf),mean(ytot(:,1:itf).^2,1),1);
+% plot(time_,fit(1)*time_+fit(2),'--k','DisplayName',['$\alpha=',num2str(fit(1)),'$']); hold on
+% plot(time_,mean(ytot.^2,1),'DisplayName','$\langle y.^2\rangle_p$');
+%
+% % r^2 displacement
+% fit = polyfit(time_(1:itf),mean(xtot(:,1:itf).^2+ytot(:,1:itf).^2,1),1);
+% plot(time_,fit(1)*time_+fit(2),'--k','DisplayName',['$\alpha=',num2str(fit(1)),'$']); hold on
+% plot(time_,mean(xtot.^2+ytot.^2,1),'DisplayName','$\langle r.^2\rangle_p$');
+% ylabel('$\langle x^2 \rangle_p$');
+% xlim(U_TIME + [0 Tfin]);
subplot(223);
for ip = 1:Np
@@ -247,14 +317,8 @@ ylabel('$y_p$');
xlim(U_TIME + [0 Tfin]);
subplot(224);
- plot(time_,mean(ytot,1)); hold on
-xlabel('time');
-ylabel('$\langle y \rangle_p$');
-xlim(U_TIME + [0 Tfin]);
+histogram(xtot(:,1),20); hold on
+histogram(xtot(:,end),20)
+xlabel('position');
+ylabel('$n$');
-if 0
- %%
- figure
-% hist(reshape(xtot,[Np*(Nt-1) 1],Np))
- hist(xtot(:,2000),Np/20)
-end
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index a6dd80cb252bfe9e60212251b99ee403c3bae014..82ec97b90edcd0b88841ac038f2288e43e2cd8fb 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -6,16 +6,16 @@ default_plots_options
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 1.8; % Density gradient drive
+K_N = 1.7; % Density gradient drive
K_T = 0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
NX = 64; % real space x-gridpoints
NY = 1; % '' y-gridpoints
-LX = 100; % Size of the squared frequency domain
+LX = 120; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
Q0 = 1.0; % safety factor
@@ -23,8 +23,8 @@ SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
-DT = 1e-2c; % Time step
+TMAX = 300; % Maximal time unit
+DT = 2e-2; % 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
@@ -32,13 +32,13 @@ SPS5D = 1; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-SIMID = 'Linear_Hallenbert'; % Name of the simulation
+SIMID = 'Linear_scan'; % Name of the simulation
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'SG';
-GKCO = 0; % gyrokinetic operator
+GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
@@ -54,9 +54,9 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
+HD_CO = 0.0; % Hyper diffusivity cutoff ratio
kmax = NX*pi/LX;% Highest fourier mode
-MU = 0.05; % Hyperdiffusivity coefficient
+MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
@@ -71,17 +71,21 @@ CURVB = 1.0;
%% PARAMETER SCANS
if 1
-%% Parameter scan over PJ
-% PA = [20];
-% JA = [10];
-PA = [4 6 8];
-JA = [2 3 4];
-DTA= DT*ones(size(JA));%./sqrt(JA);
-% DTA= DT;
-mup_ = MU_P;
-muj_ = MU_J;
+% Parameter scan over PJ
+PA = [4 6 8 10];
+JA = [2 3 4 5];
Nparam = numel(PA);
-param_name = 'PJ';
+% Parameter scan over KN
+% PA = [4]; JA = [2];
+% PMAXE = PA(1); PMAXI = PA(1);
+% JMAXE = JA(1); JMAXI = JA(1);
+% KNA = 1.5:0.05:2.5;
+% ETA = 0.25;
+% Nparam = numel(KNA);
+%
+DTA= DT*ones(1,Nparam)./sqrt(JA);
+% DTA= DT;
+param_name = 'KN';
gamma_Ni00 = zeros(Nparam,floor(NX/2)+1);
gamma_Nipj = zeros(Nparam,floor(NX/2)+1);
gamma_phi = zeros(Nparam,floor(NX/2)+1);
@@ -89,12 +93,15 @@ for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
JMAXE = JA(i); JMAXI = JA(i);
+% K_N = KNA(i); K_T = ETA*K_N;
DT = DTA(i);
setup
system(['rm fort*.90']);
% Run linear simulation
if RUN
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0; cd ../../../wk'])
+% disp([param_name,'=',num2str(K_N)]);
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0 > out.txt; cd ../../../wk']);
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
@@ -127,7 +134,7 @@ for i = 1:Nparam
end
if 1
-%% Plot
+%% Plot for PJ scan
SCALE = 1;%sqrt(2);
fig = figure; FIGNAME = 'linear_study';
plt = @(x) x;
@@ -151,6 +158,19 @@ saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']);
saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']);
end
end
+
+if 0
+%% Plot for KN scan
+fig = figure; FIGNAME = 'linear_study';
+[Y,X] = meshgrid(kx,KNA);
+pclr = pcolor(Y,X,gamma_phi);set(pclr, 'edgecolor','none'); colorbar;
+shading interp
+% imagesc(kx,KNA,gamma_phi);
+set(gca,'YDir','normal')
+colormap(bluewhitered); xlim([kx(2) kx(end)]);
+title(['$\gamma^p$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
+xlabel('$k_x$'); ylabel('$\kappa_N$');
+end
if 0
%% Space time
[YT,XT] = meshgrid(Ts3D,kx);
diff --git a/wk/load_f_gene.m b/wk/load_f_gene.m
new file mode 100644
index 0000000000000000000000000000000000000000..d31640e699109b7eb791c55e3b9f2336c1420cd8
--- /dev/null
+++ b/wk/load_f_gene.m
@@ -0,0 +1,99 @@
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_36x20/';
+folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/';
+% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/HP_fig_2b_mu_1e-1/';
+% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
+% folder = '/misc/gene_results/HP_fig_2c_mu_1e-2_muv_1e-1/';
+% folder = '/misc/gene_results/HP_fig_2c_gyroLES/';
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_128x32x48x32/';
+% coordData=loadCoord(folder,-1,1,0,'ions');
+% info=dat_file_parser(folder,-1);
+%
+% chkpt_file=[folder,'vsp.dat'];
+% fid=fopen(chkpt_file);
+%
+% [~,~,gene_data]=GetFortranStep(fid,'DOUBLE',8,'LITTLE',1);
+% fclose(fid);
+%
+% nx = info.box.nx0;
+% ny = info.box.nky0;
+% nz = info.box.nz0;
+% nv = info.box.nv0;
+% nw = info.box.nw0;
+% ns = numel(info.species);
+% Lw = info.box.lw;
+% Lv = info.box.lv;
+% s = linspace(-Lv,Lv,nv);
+% x = linspace(0,Lw,nw);
+% [XX,SS] = meshgrid(x,s);
+%
+% gridsize=(nv*nw*nz);
+% Gdata =reshape(gene_data((4*gridsize+1):(5*gridsize)),[nz nv nw]);
+
+file = 'coord.dat.h5';
+vp = h5read([folder,file],'/coord/vp');
+mu = h5read([folder,file],'/coord/mu');
+[XX,SS] = meshgrid(mu,vp);
+
+[~,iv0] = min(abs(vp));
+[~,im0] = min(abs(mu));
+
+file = 'vsp.dat.h5';
+time = h5read([folder,file],'/vsp/time');
+
+TIMES = 500:1000;
+
+fig = figure;
+
+G_t = [];
+Gdata = 0;
+for T = TIMES
+[~, it] = min(abs(time-T));
+tmp = h5read([folder,file],['/vsp/<f_>/',sprintf('%10.10d',it-1)]);
+Gdata = Gdata + tmp;
+G_t = [G_t time(it)];
+end
+Gdata = Gdata ./ numel(TIMES);
+
+if 0
+ FFe = squeeze(Gdata(1,:,:,2));
+ FFe = abs(FFe)./max(max(abs(FFe)));
+subplot(1,2,1)
+ plot(vp,FFe(:,im0)); hold on;
+subplot(1,2,2)
+ plot(mu,FFe(iv0,:)); hold on;
+
+ FFi = squeeze(Gdata(1,:,:,1));
+ FFi = abs(FFi)./max(max(abs(FFi)));
+subplot(1,2,1)
+ plot(vp,FFi(:,im0)); hold on;
+ legend('e','i')
+ xlabel('$v_\parallel, (\mu=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+
+subplot(1,2,2)
+ plot(mu,FFi(iv0,:)); hold on;
+ legend('e','i')
+ xlabel('$\mu, (v_\parallel=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+
+else
+subplot(1,2,1)
+ FFe = squeeze(Gdata(1,:,:,1));
+ FFe = abs(FFe)./max(max(abs(FFe)));
+ contour(SS,XX,FFe,128);
+% pclr = pcolor(SS,XX,FFe); set(pclr, 'edgecolor','none'); shading interp
+subplot(1,2,2)
+ FFi = squeeze(Gdata(1,:,:,2));
+ FFi = abs(FFi)./max(max(abs(FFi)));
+ contour(SS,XX,FFi,128);
+% pclr = pcolor(SS,XX,FFi); set(pclr, 'edgecolor','none'); shading interp
+end
+
+subplot(1,2,1)
+if numel(TIMES) == 1
+ title(['Gene, $t = ',num2str(time(it)),'$']);
+else
+ title([' Gene, average $t\in$[',num2str(G_t(1)),',',num2str(G_t(end)),']']);
+end
+
+clear FFi FFe tmp Gdata
+
\ No newline at end of file
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 841237cb613151baef69e6fb9b1e29cc66dd7f32..5e3214060d6566251889967215cc6c2ff6b5177e 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,10 +1,10 @@
clear all;
addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
-CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain)
+CHAIN = 0; % To chain jobs (CHAIN = n will launch n jobs in chain)
% EXECNAME = 'helaz3_dbg';
- EXECNAME = 'helaz3.03';
- SIMID = 'simulation_A_new';
+ EXECNAME = 'helaz3';
+ SIMID = 'HP_fig2b_conv';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -13,15 +13,16 @@ CLUSTER.PART = 'prod'; % dbg or prod
% CLUSTER.PART = 'dbg';
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
-CLUSTER.MEM = '128GB'; % Memory
-CLUSTER.JNAME = 'HeLaZ';% Job name
-NP_P = 2; % MPI processes along p
-NP_KX = 24; % MPI processes along kx
+CLUSTER.MEM = '16GB'; % Memory
+CLUSTER.JNAME = 'nu0_b_conv';% Job name
+NP_P = 1; % MPI processes along p
+NP_KX = 12; % MPI processes along kx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
-K_N = 1.0/0.6; % Density gradient drive (R/Ln)
-NU_HYP = 0.0;
+NU = 0.0; % Collision frequency
+K_N = 2.5; % Density gradient drive (R/Ln)
+K_T = K_N/4; % Temperature gradient
+MU = 0.1;
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
NX = 300; % Realspace x-gridpoints
@@ -35,17 +36,17 @@ EPS = 0.0; % inverse aspect ratio
P = 4;
J = 2;
%% TIME PARAMETERS
-TMAX = 10000; % Maximal time unit
-DT = 7e-3; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1; % Sampling per time unit for 2D arrays
+TMAX = 1000; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1/2; % Sampling per time unit for profiler
+SPS2D = 1/2; % Sampling per time unit for 2D arrays
SPS3D = 1/2; % Sampling per time unit for 3D arrays
SPS5D = 1/50; % Sampling per time unit for 5D arrays
-JOB2LOAD= 0; % start from t=0 if <0, else restart from outputs_$job2load
+JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load
%% OPTIONS AND NAMING
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'LR';
+CO = 'DG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
CLOS = 0; % Closure model (0: =0 truncation)
@@ -53,9 +54,10 @@ NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >
LINEARITY = 'nonlinear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
+INIT_OPT = 'phi';
+ACT_ON_MODES = 'donothing';
%% OUTPUTS
-W_DOUBLE = 1;
+W_DOUBLE = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 1; W_TEMP = 1;
@@ -76,16 +78,12 @@ KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KXEQ0 = 0; % put kx = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = NX*pi/LX;% Highest fourier mode
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-% kmaxcut = 2.5;
-MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
+MU_X = 1e-1; % Hyperdiffusivity coefficient
+MU_Y = 1e-1; % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
BCKGD0 = 0.0; % Init background
TAU = 1.0; % e/i temperature ratio
-K_T = 0.0; % Temperature gradient
K_E = 0.0; % ES '''
-INIT_PHI= 1; % Start simulation with a noisy phi and moments
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
% Compute processes distribution
diff --git a/wk/mode_growth_meter.m b/wk/mode_growth_meter.m
index 3b90ca720d00b0c61e4217d0268a39899b29f78e..5ff7721e8ae9e8fc7a1d7662b19c608d5715399a 100644
--- a/wk/mode_growth_meter.m
+++ b/wk/mode_growth_meter.m
@@ -1,42 +1,45 @@
-MODES_SELECTOR = 2; %(1:Zonal, 2: NZonal, 3: ky=kx)
-NORMALIZED = 0;
-options.K2PLOT = 0.01:0.01:1.0;
-options.TIME =4000:1:4200;
-DATA = data;
-OPTIONS = options;
+function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
+NORMALIZED = OPTIONS.NORMALIZED;
+Nma = OPTIONS.NMA; %Number moving average
t = OPTIONS.TIME;
iz = 1;
[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
FRAMES = zeros(size(OPTIONS.TIME));
+
for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
end
+FIGURE.fig = figure; set(gcf, 'Position', [100 100 1200 700])
+FIGURE.FIGNAME = 'mode_growth_meter';
+for i = 1:2
+MODES_SELECTOR = i; %(1:Zonal, 2: NZonal, 3: ky=kx)
+
if MODES_SELECTOR == 1
if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma);
else
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma);
end
kstr = 'k_x';
k = DATA.kx;
MODESTR = 'Zonal modes';
elseif MODES_SELECTOR == 2
if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma);
else
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma);
end
kstr = 'k_y';
k = DATA.ky;
MODESTR = 'NZ modes';
elseif MODES_SELECTOR == 3
if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),Nma);
else
- plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),1);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),Nma);
end
kstr = 'k_y=k_x';
k = DATA.ky;
@@ -63,16 +66,16 @@ for ik = MODES
end
%plot
-figure; set(gcf, 'Position', [100 100 1200 350])
-subplot(1,3,1)
+subplot(2,3,1+3*(i-1))
[YY,XX] = meshgrid(t,fftshift(k,1));
pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
+% pclr = imagesc(t,fftshift(k,1),abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));
% xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
title(MODESTR)
-subplot(1,3,2)
- mod2plot = round(linspace(2,numel(MODES),15));
+subplot(2,3,2+3*(i-1))
+ mod2plot = [2:OPTIONS.NMODES+1];
for i_ = mod2plot
semilogy(t,plt(DATA.PHI,MODES(i_))); hold on;
% semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
@@ -84,11 +87,13 @@ subplot(1,3,2)
xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
title('Measure time window')
-subplot(1,3,3)
+subplot(2,3,3+3*(i-1))
plot(k(MODES),gamma); hold on;
plot(k(MODES(mod2plot)),gamma(mod2plot),'x')
if MODES_SELECTOR == 1
plot(k(ikzf),gamma(ikzf),'ok');
end
xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
- title('Growth rates')
\ No newline at end of file
+ title('Growth rates')
+end
+end
\ No newline at end of file
diff --git a/wk/quick_gene_load.m b/wk/quick_gene_load.m
index 40cb7617192efc0128a863a2291c34ee8c9916cf..7fd7a569fbc8addf16b6e6c721d5b715dcee35c9 100644
--- a/wk/quick_gene_load.m
+++ b/wk/quick_gene_load.m
@@ -29,7 +29,14 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.235_64x32.txt';
% fname ='GENE_LIN_Kn_1.7_KT_0.425_nuSG_0.235_64x32.txt';
-fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.475_nu_0_mu_5e-2.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_2.0_KT_0.5_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
+fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
data_ = load([path,fname]);
figure