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 0000000..904bc6d
--- /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 1957f42..0e81e20 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,145 +1,120 @@
-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;
Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
[~,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
% Compute Gamma from ifft matlab
Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
for it = 1:numel(DATA.Ts3D)
for iz = 1:DATA.Nz
Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
.*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
end
Gx(:,:,it) = Gx(:,:,it)/DATA.Nz;
end
Gx_t_mtlb = squeeze(mean(mean(Gx,1),2));
% Compute Heat flux from ifft matlab
Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
for it = 1:numel(DATA.Ts3D)
for iz = 1:DATA.Nz
Qx(:,:,it) = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
.*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
end
Qx(:,:,it) = Qx(:,:,it)/DATA.Nz;
end
Qx_t_mtlb = squeeze(mean(mean(Qx,1),2));
% zonal vs nonzonal energies for phi(t)
E_Zmode_SK = zeros(1,Ns3D);
E_NZmode_SK = zeros(1,Ns3D);
for it = 1:numel(DATA.Ts3D)
E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(DATA.PHI(ikzf,1,1,it))^2);
E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2.*(KY~=0)))));
end
[~,its3D] = min(abs(DATA.Ts3D-tstart));
[~,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
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(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
ylabel('$\Gamma_x$')
ylim([0,5*abs(Gx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
yyaxis right
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
ylabel('$Q_x$')
ylim([0,5*abs(Qx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
%% 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
%% Zonal vs NZonal energies
subplot(312)
it0 = 1; itend = Ns3D;
trange = it0:itend;
plt1 = @(x) x;%-x(1);
plt2 = @(x) x./max((x(its3D:ite3D)));
toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
% plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
plot(plt1(DATA.Ts3D(trange)),plt2(toplot(trange)),'DisplayName','Sum $A^2$');
ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
\ No newline at end of file
diff --git a/matlab/plot/real_plot_1D.m b/matlab/plot/real_plot_1D.m
index ec3217f..d0e93cd 100644
--- a/matlab/plot/real_plot_1D.m
+++ b/matlab/plot/real_plot_1D.m
@@ -1,171 +1,175 @@
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;
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.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$';
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 = compy(toplot.FIELD(:,:,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,3,2)
+subplot(1,nplots,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 = compx(toplot.FIELD(:,:,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)]);
%% Z plot
if data.Nz > 1
options.PLAN = 'yz';
options.COMP = options.COMPX;
options.POLARPLOT = 0;
options.INTERP = 0;
toplot = process_field(data,options);
t = data.Ts3D; frames = toplot.FRAMES;
-end
switch options.COMPZ
case 'avg'
compz = @(x) mean(x,1);
ynamez = ['$\langle ',yname,'\rangle_y$'];
otherwise
compz = @(x) x(1,:);
ynamez = ['$',yname,'(z,y=0)$'];
end
-subplot(1,3,3)
+subplot(1,nplots,3)
X = data.z;
xname = '$z$';
switch options.COMPT
case 'avg'
Y = compz(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(ynamez)
xlim([min(X),max(X)]);
end
-
+end
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index 926eef8..0723b80 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -1,111 +1,110 @@
function [ FIGURE ] = show_geometry(DATA,OPTIONS)
% filtering Z pinch and torus
if DATA.Nz > 1 % Torus flux tube geometry
Nptor = DATA.Nz; Tgeom = 1;
Nturns = 1;
a = DATA.EPS; % Torus minor radius
R = 1.; % Torus major radius
q = (DATA.Nz>1)*DATA.Q0; % Flux tube safety factor
theta = linspace(-pi, pi, Nptor) ; % Poloidal angle
phi = linspace(0., 2.*pi, Nptor) ; % Toroidal angle
[t, p] = meshgrid(phi, theta);
x_tor = (R + a.*cos(p)) .* cos(t);
y_tor = (R + a.*cos(p)) .* sin(t);
z_tor = a.*sin(p);
DIMENSIONS = [50 50 1200 600];
else % Zpinch geometry
Nturns = 0.1; Tgeom = 0;
a = 0.7; % Torus minor radius
R = 0; % Torus major radius
q = 0;
theta = linspace(-Nturns*pi, Nturns*pi, 100); % Poloidal angle
phi = linspace(-0.1, 0.1, 5); % cylinder height
[t, p] = meshgrid(phi, theta);
x_tor = a.*cos(p);
z_tor = a.*sin(p);
y_tor = t;
DIMENSIONS = [50 50 numel(OPTIONS.TIME)*400 500];
end
% field line coordinates
Xfl = @(z) (R+a*cos(z)).*cos(q*z);
Yfl = @(z) (R+a*cos(z)).*sin(q*z);
Zfl = @(z) a*sin(z);
Rvec= @(z) [Xfl(z); Yfl(z); Zfl(z)];
% xvec
xX = @(z) (Xfl(z)-R*cos(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
xY = @(z) (Yfl(z)-R*sin(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
xvec= @(z) [xX(z); xY(z); xZ(z)];
% bvec
bX = @(z) Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z))./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2);
bY = @(z) (a*cos(z).*sin(q*z) + q*Xfl(z))./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2);
bZ = @(z) a*cos(z)./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2);
bvec= @(z) [bX(z); bY(z); bZ(z)];
% yvec = b times x
yX = @(z) bY(z).*xZ(z)-bZ(z).*xY(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2);
yY = @(z) bZ(z).*xX(z)-bX(z).*xZ(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2);
yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2);
yvec= @(z) [yX(z); yY(z); yZ(z)];
% Plot high res field line
% Planes plot
theta = DATA.z ; % Poloidal angle
% Plot x,y,bvec at these points
scale = 0.10;
%% Plot plane result
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)
subplot(1,numel(OPTIONS.TIME),it_)
%plot magnetic geometry
if OPTIONS.PLT_MTOPO
magnetic_topo=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
set(magnetic_topo,'edgecolor',[1 1 1]*0.7,'facecolor','none')
end
%plot field line
theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta)); hold on;
%plot vector basis
theta = DATA.z ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r');
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
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
for iz = OPTIONS.PLANES
z_ = DATA.z(iz);
Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
s=surface(Xp,Yp,Zp,FIELDS(:,:,iz)/max_);
s.EdgeColor = 'none';
colormap(bluewhitered);
% caxis([-1,1]);
end
if DATA.Nz == 1
xlim([0.65 0.75]);
ylim([-0.1 0.1]);
zlim([-0.2 0.2]);
end
end
%%
axis equal
view([1,-2,1])
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 5258ef7..5709653 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -1,451 +1,478 @@
function [ TOPLOT ] = process_field( DATA, OPTIONS )
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
%% Setup time
%%
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;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
[X,Y] = meshgrid(DATA.x,DATA.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.x);
REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.y);
REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
[X,Y] = meshgrid(DATA.kx,DATA.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
REALP = 0; COMPDIM = 2; 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_]);
Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2;
Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2;
DIMENSIONS = [100, 100, 400*Rx, 400*Ry];
ASPECT = [Rx, Ry, 1];
% Polar grid
POLARNAME = [];
if POLARPLOT
POLARNAME = 'polar';
X__ = (X+DATA.a).*cos(Y);
Y__ = (X+DATA.a).*sin(Y);
X = X__;
Y = Y__;
XNAME='X';
YNAME='Z';
DIMENSIONS = [100, 100, 500, 500];
ASPECT = [1,1,1];
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
else
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
COMPNAME = ['avg',directions{COMPDIM}];
FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}];
case 'max'
compr = @(x) max(x,[],COMPDIM);
COMPNAME = ['max',directions{COMPDIM}];
FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME];
otherwise
switch COMPDIM
case 1
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
COMPNAME = sprintf(['x=','%2.1f'],DATA.x(i));
else
COMPNAME = sprintf(['k_x=','%2.1f'],DATA.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
COMPNAME = sprintf(['y=','%2.1f'],DATA.y(i));
else
COMPNAME = sprintf(['k_y=','%2.1f'],DATA.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
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
INTERP = 0;
switch COMPDIM
case 1
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,1);
shift_y = @(x) fftshift(x,1);
case 2
process = @(x) abs((x));
shift_x = @(x) (x);
shift_y = @(x) (x);
case 3
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,2);
shift_y = @(x) fftshift(x,2);
end
end
switch OPTIONS.NAME
case '\phi'
NAME = 'phi';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'N_i^{00}'
NAME = 'Ni00';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.Ni00(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'n_e'
NAME = 'ne';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
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
tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'n_e^{NZ}'
NAME = 'neNZ';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'n_i'
NAME = 'ni';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'T_i'
NAME = 'Ti';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.TEMP_I(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case 'n_i^{NZ}'
NAME = 'niNZ';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
case '\phi^{NZ}'
NAME = 'phiNZ';
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
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)
tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
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
for it = 1:numel(FRAMES)
FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
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)))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
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
for it = 1:numel(FRAMES)
FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
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)))));
end
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 = 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
for it = 1:numel(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 = 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
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));
+ [~, 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
for it = 1:numel(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);
+ 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'
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 = 1:numel(FRAMES)
OPTIONS.T = DATA.Ts5D(FRAMES(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 = 1:numel(FRAMES)
OPTIONS.T = DATA.Ts5D(FRAMES(it));
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
end
TOPLOT.FIELD = FIELD;
TOPLOT.FRAMES = FRAMES;
TOPLOT.X = shift_x(X);
TOPLOT.Y = shift_y(Y);
TOPLOT.FIELDNAME = FIELDNAME;
TOPLOT.XNAME = XNAME;
TOPLOT.YNAME = YNAME;
TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME];
TOPLOT.DIMENSIONS= DIMENSIONS;
TOPLOT.ASPECT = ASPECT;
TOPLOT.FRAMES = FRAMES;
TOPLOT.INTERP = INTERP;
end
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 50b9e5c..36585c3 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -1,525 +1,525 @@
MODULE grid
! Grid module for spatial discretization
USE prec_const
USE basic
IMPLICIT NONE
PRIVATE
! GRID Namelist
INTEGER, PUBLIC, PROTECTED :: pmaxe = 1 ! The maximal electron Hermite-moment computed
INTEGER, PUBLIC, PROTECTED :: jmaxe = 1 ! The maximal electron Laguerre-moment computed
INTEGER, PUBLIC, PROTECTED :: pmaxi = 1 ! The maximal ion Hermite-moment computed
INTEGER, PUBLIC, PROTECTED :: jmaxi = 1 ! The maximal ion Laguerre-moment computed
INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal Laguerre-moment
INTEGER, PUBLIC, PROTECTED :: dmaxe = 1 ! The maximal full GF set of e-moments v^dmax
INTEGER, PUBLIC, PROTECTED :: dmaxi = 1 ! The maximal full GF set of i-moments v^dmax
INTEGER, PUBLIC, PROTECTED :: Nx = 16 ! Number of total internal grid points in x
REAL(dp), PUBLIC, PROTECTED :: Lx = 1._dp ! horizontal length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Ny = 16 ! Number of total internal grid points in y
REAL(dp), PUBLIC, PROTECTED :: Ly = 1._dp ! vertical length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nz = 1 ! Number of total perpendicular planes
REAL(dp), PUBLIC, PROTECTED :: Npol = 1._dp ! number of poloidal turns
INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes
INTEGER, PUBLIC, PROTECTED :: Nkx = 8 ! Number of total internal grid points in kx
REAL(dp), PUBLIC, PROTECTED :: Lkx = 1._dp ! horizontal length of the fourier box
INTEGER, PUBLIC, PROTECTED :: Nky = 16 ! Number of total internal grid points in ky
REAL(dp), PUBLIC, PROTECTED :: Lky = 1._dp ! vertical length of the fourier box
REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component
! For Orszag filter
REAL(dp), PUBLIC, PROTECTED :: two_third_kxmax
REAL(dp), PUBLIC, PROTECTED :: two_third_kymax
REAL(dp), PUBLIC :: two_third_kpmax
! 1D Antialiasing arrays (2/3 rule)
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_x
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_y
! Grids containing position in physical space
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: xarray
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: yarray
! Local and global z grids, 2D since it has to store odd and even grids
REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: zarray
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray_full
! local z weights for computing simpson rule
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: zweights_SR
REAL(dp), PUBLIC, PROTECTED :: deltax, deltay, deltaz, inv_deltaz, diff_dz_coeff
INTEGER, PUBLIC, PROTECTED :: ixs, ixe, iys, iye, izs, ize
INTEGER, PUBLIC, PROTECTED :: izgs, izge ! ghosts
LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag
INTEGER, PUBLIC :: ir,iz ! counters
! Data about parallel distribution for kx
integer(C_INTPTR_T), PUBLIC :: local_nkx, local_nky
integer(C_INTPTR_T), PUBLIC :: local_nkx_offset, local_nky_offset
INTEGER, PUBLIC :: local_nkp
! "" for p
INTEGER, PUBLIC :: local_np_e, local_np_i
INTEGER, PUBLIC :: total_np_e, total_np_i
integer(C_INTPTR_T), PUBLIC :: local_np_e_offset, local_np_i_offset
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_np_e, counts_np_i
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_np_e, displs_np_i
! "" for z
INTEGER, PUBLIC :: local_nz
INTEGER, PUBLIC :: total_nz
integer(C_INTPTR_T), PUBLIC :: local_nz_offset
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nz
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nz
! "" for j (not parallelized)
INTEGER, PUBLIC :: local_nj_e, local_nj_i
! Grids containing position in fourier space
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kxarray, kxarray_full
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, 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
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i, parray_i_full
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e, jarray_e_full
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i, jarray_i_full
INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e ! Start and end indices for pol. deg.
INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i
INTEGER, PUBLIC, PROTECTED :: ipgs_e,ipge_e, ijgs_e,ijge_e ! Ghosts start and end indices
INTEGER, PUBLIC, PROTECTED :: ipgs_i,ipge_i, ijgs_i,ijge_i
INTEGER, PUBLIC, PROTECTED :: deltape, ip0_e, ip1_e, ip2_e ! Pgrid spacing and moment 0,1,2 index
INTEGER, PUBLIC, PROTECTED :: deltapi, ip0_i, ip1_i, ip2_i
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip0_e, CONTAINS_ip0_i
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip1_e, CONTAINS_ip1_i
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip2_e, CONTAINS_ip2_i
! Usefull inverse numbers
REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz
! Public Functions
PUBLIC :: init_1Dgrid_distr
PUBLIC :: set_pgrid, set_jgrid
PUBLIC :: set_kxgrid, set_kygrid, set_zgrid
PUBLIC :: grid_readinputs, grid_outputinputs
PUBLIC :: bare, bari
! Precomputations
real(dp), PUBLIC, PROTECTED :: pmaxe_dp, pmaxi_dp, jmaxe_dp,jmaxi_dp
CONTAINS
SUBROUTINE grid_readinputs
! Read the input parameters
USE prec_const
IMPLICIT NONE
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
Nx, Lx, Ny, Ly, Nz, Npol, SG
READ(lu_in,grid)
!! Compute the maximal degree of full GF moments set
! i.e. : all moments N_a^pj s.t. p+2j<=d are simulated (see GF closure)
dmaxe = min(pmaxe,2*jmaxe+1)
dmaxi = min(pmaxi,2*jmaxi+1)
! If no parallel dim (Nz=1), the moment hierarchy is separable between odds and even P
!! and since the energy is injected in P=0 and P=2 for density/temperature gradients
!! there is no need of simulating the odd p which will only be damped.
!! We define in this case a grid Parray = 0,2,4,...,Pmax i.e. deltap = 2 instead of 1
!! to spare computation
IF(Nz .EQ. 1) THEN
deltape = 2; deltapi = 2;
pp2 = 1; ! index p+2 is ip+1
SG = .FALSE.
ELSE
deltape = 1; deltapi = 1;
pp2 = 2; ! index p+2 is ip+1
ENDIF
! Usefull precomputations
inv_Nx = 1._dp/REAL(Nx,dp)
inv_Ny = 1._dp/REAL(Ny,dp)
END SUBROUTINE grid_readinputs
SUBROUTINE init_1Dgrid_distr
! write(*,*) Nx
local_nky = (Ny/2+1)/num_procs_ky
! write(*,*) local_nkx
local_nky_offset = rank_ky*local_nky
if (rank_ky .EQ. num_procs_ky-1) local_nky = (Ny/2+1)-local_nky_offset
END SUBROUTINE init_1Dgrid_distr
SUBROUTINE set_pgrid
USE prec_const
IMPLICIT NONE
INTEGER :: ip, istart, iend, in
! Total number of Hermite polynomials we will evolve
total_np_e = (Pmaxe/deltape) + 1
total_np_i = (Pmaxi/deltapi) + 1
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(parray_e_full(1:total_np_e))
ALLOCATE(parray_i_full(1:total_np_i))
! P
DO ip = 1,total_np_e; parray_e_full(ip) = (ip-1)*deltape; END DO
DO ip = 1,total_np_i; parray_i_full(ip) = (ip-1)*deltapi; END DO
!! Parallel data distribution
! Local data distribution
CALL decomp1D(total_np_e, num_procs_p, rank_p, ips_e, ipe_e)
CALL decomp1D(total_np_i, num_procs_p, rank_p, ips_i, ipe_i)
local_np_e = ipe_e - ips_e + 1
local_np_i = ipe_i - ips_i + 1
! Ghosts boundaries
ipgs_e = ips_e - 2/deltape; ipge_e = ipe_e + 2/deltape;
ipgs_i = ips_i - 2/deltapi; ipge_i = ipe_i + 2/deltapi;
! List of shift and local numbers between the different processes (used in scatterv and gatherv)
ALLOCATE(counts_np_e (1:num_procs_p))
ALLOCATE(counts_np_i (1:num_procs_p))
ALLOCATE(displs_np_e (1:num_procs_p))
ALLOCATE(displs_np_i (1:num_procs_p))
DO in = 0,num_procs_p-1
CALL decomp1D(total_np_e, num_procs_p, in, istart, iend)
counts_np_e(in+1) = iend-istart+1
displs_np_e(in+1) = istart-1
CALL decomp1D(total_np_i, num_procs_p, in, istart, iend)
counts_np_i(in+1) = iend-istart+1
displs_np_i(in+1) = istart-1
ENDDO
! local grid computation
CONTAINS_ip0_e = .FALSE.
CONTAINS_ip1_e = .FALSE.
CONTAINS_ip2_e = .FALSE.
CONTAINS_ip0_i = .FALSE.
CONTAINS_ip1_i = .FALSE.
CONTAINS_ip2_i = .FALSE.
ALLOCATE(parray_e(ipgs_e:ipge_e))
ALLOCATE(parray_i(ipgs_i:ipge_i))
DO ip = ipgs_e,ipge_e
parray_e(ip) = (ip-1)*deltape
! Storing indices of particular degrees for fluid moments computations
IF(parray_e(ip) .EQ. 0) THEN
ip0_e = ip
CONTAINS_ip0_e = .TRUE.
ENDIF
IF(parray_e(ip) .EQ. 1) THEN
ip1_e = ip
CONTAINS_ip1_e = .TRUE.
ENDIF
IF(parray_e(ip) .EQ. 2) THEN
ip2_e = ip
CONTAINS_ip2_e = .TRUE.
ENDIF
END DO
DO ip = ipgs_i,ipge_i
parray_i(ip) = (ip-1)*deltapi
! Storing indices of particular degrees for fluid moments computations
IF(parray_i(ip) .EQ. 0) THEN
ip0_i = ip
CONTAINS_ip0_i = .TRUE.
ENDIF
IF(parray_i(ip) .EQ. 1) THEN
ip1_i = ip
CONTAINS_ip1_i = .TRUE.
ENDIF
IF(parray_i(ip) .EQ. 2) THEN
ip2_i = ip
CONTAINS_ip2_i = .TRUE.
ENDIF
END DO
!DGGK operator uses moments at index p=2 (ip=3) for the p=0 term so the
! process that contains ip=1 MUST contain ip=3 as well for both e and i.
IF(((ips_e .EQ. ip0_e) .OR. (ips_i .EQ. ip0_e)) .AND. ((ipe_e .LT. ip2_e) .OR. (ipe_i .LT. ip2_i)))&
WRITE(*,*) "Warning : distribution along p may not work with DGGK"
! Precomputations
pmaxe_dp = real(pmaxe,dp)
pmaxi_dp = real(pmaxi,dp)
END SUBROUTINE set_pgrid
SUBROUTINE set_jgrid
USE prec_const
IMPLICIT NONE
INTEGER :: ij
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(jarray_e_full(1:jmaxe+1))
ALLOCATE(jarray_i_full(1:jmaxi+1))
! J
DO ij = 1,jmaxe+1; jarray_e_full(ij) = (ij-1); END DO
DO ij = 1,jmaxi+1; jarray_i_full(ij) = (ij-1); END DO
! Local data
ijs_e = 1; ije_e = jmaxe + 1
ijs_i = 1; ije_i = jmaxi + 1
! Ghosts boundaries
ijgs_e = ijs_e - 1; ijge_e = ije_e + 1;
ijgs_i = ijs_i - 1; ijge_i = ije_i + 1;
! Local number of J
local_nj_e = ijge_e - ijgs_e + 1
local_nj_i = ijge_i - ijgs_i + 1
ALLOCATE(jarray_e(ijgs_e:ijge_e))
ALLOCATE(jarray_i(ijgs_i:ijge_i))
DO ij = ijgs_e,ijge_e; jarray_e(ij) = ij-1; END DO
DO ij = ijgs_i,ijge_i; jarray_i(ij) = ij-1; END DO
! Precomputations
maxj = MAX(jmaxi, jmaxe)
jmaxe_dp = real(jmaxe,dp)
jmaxi_dp = real(jmaxi,dp)
END SUBROUTINE set_jgrid
SUBROUTINE set_kygrid
USE prec_const
USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_
Nky = Ny/2+1 ! Defined only on positive kx since fields are real
! Grid spacings
IF (Ny .EQ. 1) THEN
deltaky = 0._dp
ky_max = 0._dp
ky_min = 0._dp
ELSE
deltaky = 2._dp*PI/Ly
ky_max = Nky*deltaky
ky_min = deltaky
ENDIF
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(kyarray_full(1:Nky))
DO iky = 1,Nky
kyarray_full(iky) = REAL(iky-1,dp) * deltaky
END DO
!! Parallel distribution
ikys = local_nky_offset + 1
ikye = ikys + local_nky - 1
ALLOCATE(kyarray(ikys:ikye))
local_kymax = 0._dp
! Creating a grid ordered as dk*(0 1 2 3)
DO iky = ikys,ikye
kyarray(iky) = REAL(iky-1,dp) * deltaky
! Finding kx=0
IF (kyarray(iky) .EQ. 0) THEN
iky_0 = iky
contains_ky0 = .true.
ENDIF
! Finding local kxmax value
IF (ABS(kyarray(iky)) .GT. local_kymax) THEN
local_kymax = ABS(kyarray(iky))
ENDIF
! Finding kxmax idx
IF (kyarray(iky) .EQ. ky_max) THEN
iky_max = iky
contains_kymax = .true.
ENDIF
END DO
! Orszag 2/3 filter
two_third_kymax = 2._dp/3._dp*deltaky*(Nky-1)
ALLOCATE(AA_y(ikys:ikye))
DO iky = ikys,ikye
IF ( (kyarray(iky) .LT. two_third_kymax) .OR. (LINEARITY .EQ. 'linear')) THEN
AA_y(iky) = 1._dp;
ELSE
AA_y(iky) = 0._dp;
ENDIF
END DO
END SUBROUTINE set_kygrid
SUBROUTINE set_kxgrid
USE prec_const
USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_, counter
Nkx = Nx;
! Local data
! Start and END indices of grid
ikxs = 1
ikxe = Nkx
local_nkx = ikxe - ikxs + 1
ALLOCATE(kxarray(ikxs:ikxe))
ALLOCATE(kxarray_full(1:Nkx))
IF (Nx .EQ. 1) THEN ! "cancel" y dimension
deltakx = 1._dp
kxarray(1) = 0._dp
ikx_0 = 1
contains_kx0 = .true.
kx_max = 0._dp
ikx_max = 1
kx_min = 0._dp
kxarray_full(1) = 0._dp
local_kxmax = 0._dp
ELSE ! Build apprpopriate grid
deltakx = 2._dp*PI/Lx
kx_max = (Ny/2)*deltakx
kx_min = deltakx
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
local_kxmax = 0._dp
DO ikx = ikxs,ikxe
kxarray(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
if (ikx .EQ. Nx/2+1) kxarray(ikx) = -kxarray(ikx)
! Finding kx=0
IF (kxarray(ikx) .EQ. 0) THEN
ikx_0 = ikx
contains_kx0 = .true.
ENDIF
! Finding local kxmax
IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
local_kxmax = ABS(kxarray(ikx))
ENDIF
! Finding kxmax
IF (kxarray(ikx) .EQ. kx_max) ikx_max = ikx
END DO
! Build the full grids on process 0 to diagnose it without comm
! kx
DO ikx = 1,Nkx
kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
ENDIF
! Orszag 2/3 filter
two_third_kxmax = 2._dp/3._dp*deltakx*(Nkx/2-1);
ALLOCATE(AA_x(ikxs:ikxe))
DO ikx = ikxs,ikxe
IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. &
(kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
AA_x(ikx) = 1._dp;
ELSE
AA_x(ikx) = 0._dp;
ENDIF
END DO
END SUBROUTINE set_kxgrid
SUBROUTINE set_zgrid
USE prec_const
IMPLICIT NONE
INTEGER :: i_, fid
REAL :: grid_shift, Lz
INTEGER :: ip, istart, iend, in
total_nz = Nz
! Length of the flux tube (in ballooning angle)
Lz = 2_dp*pi*Npol
! Z stepping (#interval = #points since periodic)
deltaz = Lz/REAL(Nz,dp)
inv_deltaz = 1._dp/deltaz
diff_dz_coeff = (deltaz/2._dp)**2
IF (SG) THEN
grid_shift = deltaz/2._dp
ELSE
grid_shift = 0._dp
ENDIF
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(zarray_full(1:Nz))
IF (Nz .EQ. 1) Npol = 0
DO iz = 1,total_nz
zarray_full(iz) = REAL(iz-1,dp)*deltaz - PI*REAL(Npol,dp)
END DO
!! Parallel data distribution
! Local data distribution
CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
local_nz = ize - izs + 1
! Ghosts boundaries (depend on the order of z operators)
IF(Nz .EQ. 1) THEN
izgs = izs; izge = ize;
ELSEIF(Nz .GE. 4) THEN
izgs = izs - 2; izge = ize + 2;
ELSE
ERROR STOP 'Error stop: Nz is not appropriate!!'
ENDIF
! List of shift and local numbers between the different processes (used in scatterv and gatherv)
ALLOCATE(counts_nz (1:num_procs_z))
ALLOCATE(displs_nz (1:num_procs_z))
DO in = 0,num_procs_z-1
CALL decomp1D(total_nz, num_procs_z, in, istart, iend)
counts_nz(in+1) = iend-istart+1
displs_nz(in+1) = istart-1
ENDDO
! Local z array
ALLOCATE(zarray(izgs:izge,0:1))
DO iz = izgs,izge
IF(iz .EQ. 0) THEN
zarray(iz,0) = zarray_full(total_nz)
zarray(iz,1) = zarray_full(total_nz) + grid_shift
ELSEIF(iz .EQ. -1) THEN
zarray(iz,0) = zarray_full(total_nz-1)
zarray(iz,1) = zarray_full(total_nz-1) + grid_shift
ELSEIF(iz .EQ. total_nz + 1) THEN
zarray(iz,0) = zarray_full(1)
zarray(iz,1) = zarray_full(1) + grid_shift
ELSEIF(iz .EQ. total_nz + 2) THEN
zarray(iz,0) = zarray_full(2)
zarray(iz,1) = zarray_full(2) + grid_shift
ELSE
zarray(iz,0) = zarray_full(iz)
zarray(iz,1) = zarray_full(iz) + grid_shift
ENDIF
ENDDO
! Weitghs for Simpson rule
ALLOCATE(zweights_SR(izs:ize))
DO iz = izs,ize
IF((iz .EQ. 1) .OR. (iz .EQ. Nz)) THEN
zweights_SR(iz) = 1._dp
ELSEIF(MODULO(iz-1,2)) THEN
zweights_SR(iz) = 4._dp
ELSE
zweights_SR(iz) = 2._dp
ENDIF
ENDDO
END SUBROUTINE set_zgrid
SUBROUTINE grid_outputinputs(fidres, str)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach
USE prec_const
IMPLICIT NONE
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
CALL attach(fidres, TRIM(str), "pmaxe", pmaxe)
CALL attach(fidres, TRIM(str), "jmaxe", jmaxe)
CALL attach(fidres, TRIM(str), "pmaxi", pmaxi)
CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
CALL attach(fidres, TRIM(str), "Nx", Nx)
CALL attach(fidres, TRIM(str), "Lx", Lx)
CALL attach(fidres, TRIM(str), "Ny", Ny)
CALL attach(fidres, TRIM(str), "Ly", Ly)
CALL attach(fidres, TRIM(str), "Nz", Nz)
CALL attach(fidres, TRIM(str), "Nkx", Nkx)
CALL attach(fidres, TRIM(str), "Lkx", Lkx)
CALL attach(fidres, TRIM(str), "Nky", Nky)
CALL attach(fidres, TRIM(str), "Lky", Lky)
CALL attach(fidres, TRIM(str), "SG", SG)
END SUBROUTINE grid_outputinputs
FUNCTION bare(p_,j_)
IMPLICIT NONE
INTEGER :: bare, p_, j_
bare = (jmaxe+1)*p_ + j_ + 1
END FUNCTION
FUNCTION bari(p_,j_)
IMPLICIT NONE
INTEGER :: bari, p_, j_
bari = (jmaxi+1)*p_ + j_ + 1
END FUNCTION
SUBROUTINE decomp1D( n, numprocs, myid, s, e )
INTEGER :: n, numprocs, myid, s, e
INTEGER :: nlocal
INTEGER :: deficit
nlocal = n / numprocs
s = myid * nlocal + 1
deficit = MOD(n,numprocs)
s = s + MIN(myid,deficit)
IF (myid .LT. deficit) nlocal = nlocal + 1
e = s + nlocal - 1
IF (e .GT. n .OR. myid .EQ. numprocs-1) e = n
END SUBROUTINE decomp1D
END MODULE grid
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index eba8d51..2add3a7 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,192 +1,194 @@
addpath(genpath([helazdir,'matlab'])) % ... add
addpath(genpath([helazdir,'matlab/plot'])) % ... add
addpath(genpath([helazdir,'matlab/compute'])) % ... add
addpath(genpath([helazdir,'matlab/load'])) % ... add
%% Load the results
LOCALDIR = [helazdir,'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
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
disp('Plots')
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
if 0
%% statistical transport averaging
options.T = [16000 17000];
fig = statistical_transport_averaging(data,options);
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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')
end
if 0
%% 2D snapshots
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME = 'N_i^{00}';
% 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)
end
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);
save_figure(data,fig)
end
if 0
%% 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 = 5500;
options.CTR = 1;
options.ONED = 0;
fig = plot_fa(data,options);
save_figure(data,fig)
end
if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
options.P2J = 0;
options.ST = 0;
options.NORMALIZED = 0;
fig = show_moments_spectrum(data,options);
save_figure(data,fig)
end
if 0
%% Time averaged spectrum
options.TIME = 1000:5000;
options.NORM = 0;
options.NAME = '\phi';
options.NAME = 'n_i';
options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
options.COMPXY = 0;
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
% save_figure(data,fig)
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)
end
if 0
%% Mode evolution
options.NORMALIZED = 1;
options.K2PLOT = 1;
options.TIME = 5:1:15;
options.NMA = 1;
options.NMODES = 5;
options.iz = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig)
end
if 0
%% 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
%% Metric infos
fig = plot_metric(data);
end
if 0
%% linear growth rate for 3D fluxtube
trange = [0 100];
nplots = 1;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
end
if 0
%% linear growth rate for 3D Zpinch
trange = [5 15];
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
options.kzkx = 0;
options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 3970fcc..b3cf14f 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -1,12 +1,12 @@
% Directory of the code "mypathtoHeLaZ/HeLaZ/"
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 700916a..6a08d4a 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -1,86 +1,87 @@
% 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
if 0
%% 2D snapshots
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xz';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [100 300 900];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig)
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = 'v_y';
% 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
if 0
%% Geometry
names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
'$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
'$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
figure;
subplot(311)
for i = 1:6
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
subplot(312)
for i = 7:10
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
subplot(313)
for i = 11:16
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
end