From cd9c3e804ccbd07d9fe0e3a3f84d17404857dedc Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 3 May 2022 17:32:44 +0200
Subject: [PATCH] update scripts for new ky kx rep
---
.../invert_kxky_to_kykx_gene_results.m | 21 ++++
.../plot_radial_transport_and_spacetime.m | 65 ++++-------
matlab/plot/real_plot_1D.m | 26 +++--
matlab/plot/show_geometry.m | 7 +-
matlab/process_field.m | 107 +++++++++++-------
src/grid_mod.F90 | 10 +-
wk/analysis_3D.m | 34 +++---
wk/analysis_header.m | 4 +-
wk/gene_analysis_3D.m | 21 ++--
9 files changed, 162 insertions(+), 133 deletions(-)
create mode 100644 matlab/compute/invert_kxky_to_kykx_gene_results.m
diff --git a/matlab/compute/invert_kxky_to_kykx_gene_results.m b/matlab/compute/invert_kxky_to_kykx_gene_results.m
new file mode 100644
index 00000000..904bc6d0
--- /dev/null
+++ b/matlab/compute/invert_kxky_to_kykx_gene_results.m
@@ -0,0 +1,21 @@
+function [ data ] = invert_kxky_to_kykx_gene_results( data )
+%to go from nkx nky Gene representation to nky nkx HeLaZ one
+
+rotate = @(x) permute(x,[2 1 3 4]);
+
+data.PHI = rotate(data.PHI);
+data.DENS_I = rotate(data.DENS_I);
+data.TPER_I = rotate(data.TPER_I);
+data.TPAR_I = rotate(data.TPAR_I);
+data.TEMP_I = rotate(data.TEMP_I);
+
+data.Ny = data.Nky*2-1;
+data.Nx = data.Nkx;
+
+dkx = data.kx(2); dky = data.ky(2);
+Lx = 2*pi/dkx; Ly = 2*pi/dky;
+x = linspace(-Lx/2,Lx/2,data.Nx+1); x = x(1:end-1);
+y = linspace(-Ly/2,Ly/2,data.Ny+1); y = y(1:end-1);
+data.x = x; data.y = y; data.Lx = Lx; data.Ly = Ly;
+end
+
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 1957f424..0e81e20b 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,6 +1,6 @@
-function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname,Nmvm, COMPZ)
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
%Compute steady radial transport
- tend = TAVG_1; tstart = TAVG_0;
+ tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
[~,ite0D] = min(abs(DATA.Ts0D-tend));
SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2;
@@ -11,7 +11,6 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
- z = DATA.z;
%% computations
@@ -47,7 +46,7 @@ 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);
+mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
subplot(311)
yyaxis left
@@ -70,58 +69,34 @@ mvm = @(x) movmean(x,Nmvm);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
- [KY, KX] = meshgrid(DATA.ky, DATA.kx);
- plt = @(x) mean(x(:,:,1,:),1);
+ [KX, KY] = meshgrid(DATA.kx, DATA.ky);
+ plt = @(x) mean(x(:,:,:),1);
kycut = max(DATA.ky);
kxcut = max(DATA.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
- switch COMPZ
- case 'avg'
- compz = @(f) mean(f,3);
- nameL = '$\langle$ ';
- nameR = ' $\rangle_{z}$';
- otherwise
- compz = @(f) f(:,:,COMPZ);
- nameL = '';
- nameR = [' $(z=',num2str(z(COMPZ)),')$'];
- end
- switch stfname
- case 'phi'
- phi = zeros(DATA.Ny,DATA.Nx,1,Ns3D);
- for it = 1:numel(DATA.Ts3D)
- phi(:,:,1,it) = ifourier_GENE(compz(DATA.PHI(:,:,:,it)));
- end
- f2plot = phi; fname = '$\phi$';
- case 'v_y'
- dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
- for it = 1:numel(DATA.Ts3D)
- dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*compz(DATA.PHI(:,:,:,it)).*LP);
- 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.*compz(DATA.PHI(:,:,:,it)).*LP);
- 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.*compz(DATA.PHI(:,:,:,it)).*LP);
- end
- f2plot = dx2phi; fname = '$\langle \partial_x^2\phi\rangle_y$';
- end
- clim = max(max(abs(plt(f2plot(:,:,1,its3D:ite3D)))));
+
+ OPTIONS.NAME = OPTIONS.ST_FIELD;
+ OPTIONS.PLAN = 'xy';
+ OPTIONS.COMP = 'avg';
+ OPTIONS.TIME = DATA.Ts3D;
+ OPTIONS.POLARPLOT = 0;
+ toplot = process_field(DATA,OPTIONS);
+ f2plot = toplot.FIELD;
+
+ clim = max(max(max(abs(plt(f2plot(:,:,its3D:ite3D))))));
subplot(313)
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
- legend([nameL,fname,nameR]) %colorbar;
+ legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
caxis(clim*[-1 1]);
cmap = bluewhitered(256);
colormap(cmap)
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$');
- if strcmp(stfname,'Gx')
+ if OPTIONS.INTERP
+ shading interp
+ end
+ if strcmp(OPTIONS.ST_FIELD,'Gx')
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
diff --git a/matlab/plot/real_plot_1D.m b/matlab/plot/real_plot_1D.m
index ec3217ff..d0e93cd5 100644
--- a/matlab/plot/real_plot_1D.m
+++ b/matlab/plot/real_plot_1D.m
@@ -1,4 +1,9 @@
function [ FIGURE ] = real_plot_1D( data, options )
+nplots = 0;
+if data.Nx > 1; nplots = nplots + 1; end;
+if data.Ny > 1; nplots = nplots + 1; end;
+if data.Nz > 1; nplots = nplots + 1; end;
+
options.PLAN = 'xy';
options.COMP = options.COMPZ;
@@ -33,25 +38,25 @@ end
switch options.COMPX
case 'avg'
- compx = @(x) mean(x,1);
+ compx = @(x) mean(x,2);
ynamex = ['$\langle ',yname,'\rangle_x$'];
otherwise
- compx = @(x) x(1,:);
+ compx = @(x) x(:,1);
ynamex = ['$',yname,'(y,x=0)$'];
end
switch options.COMPY
case 'avg'
- compy = @(x) mean(x,2);
+ compy = @(x) mean(x,1);
ynamey = ['$\langle ',yname,'\rangle_y$'];
otherwise
- compy = @(x) x(:,1);
+ compy = @(x) x(1,:);
ynamey = ['$',yname,'(x,y=0)$'];
end
set(gcf, 'Position', [20 50 1000 500])
-subplot(1,3,1)
+subplot(1,nplots,1)
X = data.x;
xname = '$x$';
@@ -67,7 +72,7 @@ subplot(1,3,1)
legend('show')
otherwise
for it = 1:numel(toplot.FRAMES)
- Y = compy(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = compy(toplot.FIELD(:,:,it));
Y = squeeze(Y);
if options.NORM
Y = Y./max(abs(Y));
@@ -81,7 +86,7 @@ subplot(1,3,1)
xlabel(xname); ylabel(ynamey)
xlim([min(X),max(X)]);
-subplot(1,3,2)
+subplot(1,nplots,2)
X = data.y;
xname = '$y$';
@@ -97,7 +102,7 @@ subplot(1,3,2)
legend('show')
otherwise
for it = 1:numel(toplot.FRAMES)
- Y = compx(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = compx(toplot.FIELD(:,:,it));
Y = squeeze(Y);
if options.NORM
Y = Y./max(abs(Y));
@@ -123,7 +128,6 @@ if data.Nz > 1
toplot = process_field(data,options);
t = data.Ts3D; frames = toplot.FRAMES;
-end
switch options.COMPZ
case 'avg'
@@ -134,7 +138,7 @@ switch options.COMPZ
ynamez = ['$',yname,'(z,y=0)$'];
end
-subplot(1,3,3)
+subplot(1,nplots,3)
X = data.z;
xname = '$z$';
@@ -168,4 +172,4 @@ subplot(1,3,3)
xlim([min(X),max(X)]);
end
-
+end
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index 926eef86..0723b805 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -56,9 +56,9 @@ scale = 0.10;
OPTIONS.POLARPLOT = 0;
OPTIONS.PLAN = 'xy';
r_o_R = DATA.rho_o_R;
-[Y,X] = meshgrid(r_o_R*DATA.y,r_o_R*DATA.x);
+[X,Y] = meshgrid(r_o_R*DATA.x,r_o_R*DATA.y);
max_ = 0;
-FIELDS = zeros(DATA.Nx,DATA.Ny,DATA.Nz);
+FIELDS = zeros(DATA.Ny,DATA.Nx,DATA.Nz);
FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.PARAMS]; set(gcf, 'Position', DIMENSIONS)
for it_ = 1:numel(OPTIONS.TIME)
@@ -79,11 +79,10 @@ subplot(1,numel(OPTIONS.TIME),it_)
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b');
xlabel('X');ylabel('Y');zlabel('Z');
%Plot time dependent data
- [~,it] = min(abs(OPTIONS.TIME(it_)-DATA.Ts3D));
for iz = OPTIONS.PLANES
OPTIONS.COMP = iz;
toplot = process_field(DATA,OPTIONS);
- FIELDS(:,:,iz) = toplot.FIELD(:,:,it);
+ FIELDS(:,:,iz) = toplot.FIELD(:,:,it_);
tmp = max(max(abs(FIELDS(:,:,iz))));
if (tmp > max_) max_ = tmp;
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 5258ef7c..57096536 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -7,9 +7,9 @@ FRAMES = zeros(size(OPTIONS.TIME));
for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
end
-
+FRAMES = unique(FRAMES);
%% Setup the plot geometry
-[KY,~] = meshgrid(DATA.ky,DATA.kx);
+[~,KY] = meshgrid(DATA.kx,DATA.ky);
directions = {'x','y','z'};
Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
@@ -110,7 +110,7 @@ end
switch REALP
case 1 % Real space plot
INTERP = OPTIONS.INTERP;
- process = @(x) real(fftshift(ifft2(x,Ny,Nx)));
+ process = @(x) real(fftshift(ifourier_GENE(x)));
shift_x = @(x) x;
shift_y = @(x) x;
case 0 % Frequencies plot
@@ -141,9 +141,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -161,9 +161,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -181,9 +181,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -194,7 +194,7 @@ switch OPTIONS.NAME
end
case 'k^2n_e'
NAME = 'k2ne';
- [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ [KX, KY] = meshgrid(DATA.kx, DATA.ky);
if COMPDIM == 3
for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
@@ -204,9 +204,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -224,9 +224,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -244,9 +244,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -264,9 +264,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -284,9 +284,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -304,9 +304,9 @@ switch OPTIONS.NAME
end
else
if REALP
- tmp = zeros(Nx,Ny,Nz);
+ tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
@@ -317,11 +317,11 @@ switch OPTIONS.NAME
end
case 'v_y'
NAME = 'vy';
- [~, KX] = meshgrid(DATA.ky, DATA.kx);
- vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ [KX, ~] = meshgrid(DATA.kx, DATA.ky);
+ vE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Ny,DATA.Nx));
end
end
if REALP % plot in real space
@@ -330,7 +330,7 @@ switch OPTIONS.NAME
end
else % Plot spectrum
process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it)))));
@@ -340,11 +340,11 @@ switch OPTIONS.NAME
end
case 'v_x'
NAME = 'vx';
- [KY, ~] = meshgrid(DATA.ky, DATA.kx);
- vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ [~, KY] = meshgrid(DATA.kx, DATA.ky);
+ vE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
end
if REALP % plot in real space
@@ -353,7 +353,7 @@ switch OPTIONS.NAME
end
else % Plot spectrum
process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))));
@@ -367,7 +367,7 @@ switch OPTIONS.NAME
vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifourier_GENE(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
end
if REALP % plot in real space
@@ -376,8 +376,8 @@ switch OPTIONS.NAME
end
else % Plot spectrum
process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
- for it = 1:numel(FRAMES)
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
+ for it = 1:FRAMES
for iz = 1:numel(DATA.z)
tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
@@ -386,12 +386,12 @@ switch OPTIONS.NAME
end
case '\Gamma_x'
NAME = 'Gx';
- [KY, ~] = meshgrid(DATA.ky, DATA.kx);
- Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ [~, KY] = meshgrid(DATA.kx, DATA.ky);
+ Gx = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny)))...
- .*real((ifft2(DATA.DENS_I(:,:,iz,FRAMES(it)),DATA.Nx,DATA.Ny)));
+ Gx(:,:,iz,it) = real((ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))))))...
+ .*real((ifourier_GENE(DATA.DENS_I(:,:,iz,FRAMES(it)))));
end
end
if REALP % plot in real space
@@ -402,12 +402,39 @@ switch OPTIONS.NAME
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);
+ tmp = zeros(DATA.Ny,DATA.Nx,Nz);
+ for it = 1:numel(FRAMES)
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Ny,DATA.Nx)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
+ end
+ end
+ case 'Q_x'
+ NAME = 'Qx';
+ [~, KY] = meshgrid(DATA.kx, DATA.ky);
+ Qx = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
+ for it = 1:numel(FRAMES)
+ for iz = 1:DATA.Nz
+ Qx(:,:,iz,it) = real(ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))))...
+ .*real(ifourier_GENE(DATA.DENS_I(:,:,iz,FRAMES(it))))...
+ .*real(ifourier_GENE(DATA.TEMP_I(:,:,iz,FRAMES(it))));
+ end
+ end
+ if REALP % plot in real space
+ for it = 1:numel(FRAMES)
+ FIELD(:,:,it) = squeeze(compr(Qx(:,:,:,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.Ny,DATA.Nx,Nz);
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ tmp(:,:,iz) = process(squeeze(fft2(Qx(:,:,iz,it),DATA.Ny,DATA.Nx)));
end
- FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
+ FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
end
end
case 'f_i'
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 50b9e5c7..36585c36 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -74,11 +74,11 @@ MODULE grid
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
- INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2, eo ! counters
- LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index
- LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index
- LOGICAL, PUBLIC, PROTECTED :: contains_kymax = .false. ! flag if the proc contains kx=kxmax index
+ INTEGER, PUBLIC, PROTECTED :: ikx_0, iky_0, ikx_max, iky_max ! Indices of k-grid origin and max
+ INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2, eo ! counters
+ LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index
+ LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index
+ LOGICAL, PUBLIC, PROTECTED :: contains_kymax = .false. ! flag if the proc contains kx=kxmax index
! Grid containing the polynomials degrees
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e, parray_e_full
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index eba8d519..2add3a76 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -21,10 +21,12 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 0.8*data.Ts3D(end); TAVG_1 = data.Ts3D(end); % Averaging times duration
-compz = 'avg';
-% chose your field to plot in spacetime diag (uzf,szf,Gx)
-fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi',1,compz);
+options.TAVG_0 = 0.8*data.Ts3D(end);
+options.TAVG_1 = data.Ts3D(end); % Averaging times duration
+options.NMVA = 1; % Moving average for time traces
+options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.INTERP = 1;
+fig = plot_radial_transport_and_spacetime(data,options);
save_figure(data,fig)
end
@@ -39,18 +41,18 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
-% options.NAME = '\Gamma_x';
+options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = dat.Ts5D;
-options.TIME = 100:1:200;
+options.TIME = 00:1:200;
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -68,11 +70,11 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 1;
-options.TIME = [1];
+options.TIME = [150];
data.a = data.EPS * 1000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -82,8 +84,8 @@ if 0
%% 3D plot on the geometry
options.INTERP = 1;
options.NAME = 'n_i';
-options.PLANES = 1:3:15;
-options.TIME = [100];
+options.PLANES = 1;
+options.TIME = [100 200];
options.PLT_MTOPO = 1;
data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
@@ -133,16 +135,16 @@ end
if 0
%% 1D real plot
-options.TIME = 20;
+options.TIME = [50 100 200];
options.NORM = 0;
options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME ='\Gamma_x';
% options.NAME ='s_y';
-options.COMPX = 1;
-options.COMPY = 1;
+options.COMPX = 'avg';
+options.COMPY = 'avg';
options.COMPZ = 1;
-options.COMPT = 'avg';
+options.COMPT = 1;
options.MOVMT = 1;
fig = real_plot_1D(data,options);
% save_figure(data,fig)
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 3970fcc8..b3cf14ff 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -3,10 +3,10 @@ helazdir = '/home/ahoffman/HeLaZ/';
% Directory of the simulation (from results)
% if 1% Local results
outfile ='';
-outfile ='quick_run/64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
+% outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
% outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
% outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
-% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
+outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/';
% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/';
% outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1';
run analysis_3D
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index 700916aa..6a08d4a4 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -1,21 +1,22 @@
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
-folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.2/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
gene_data = load_gene_data(folder);
-gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
+% gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
+gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 500; TAVG_1 = 700; % Averaging times duration
-% chose your field to plot in spacetime diag (uzf,szf,Gx)
-field = 'phi';
-compz = 'avg';
-nmvm = 1;
-fig = plot_radial_transport_and_spacetime(gene_data,TAVG_0,TAVG_1,field,nmvm,compz);
+options.TAVG_0 = 0.8*gene_data.Ts3D(end);
+options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
+options.NMVA = 1; % Moving average for time traces
+options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
+options.INTERP = 1;
+fig = plot_radial_transport_and_spacetime(gene_data,options);
% save_figure(data,fig)
end
@@ -50,11 +51,11 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xz';
+options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = 400:700;
+options.TIME = 000:200;
gene_data.a = data.EPS * 2000;
create_film(gene_data,options,'.gif')
end
--
GitLab