From 1a73dc547b1599b058085e4a69a2363b59494434 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 29 Apr 2022 15:13:30 +0200
Subject: [PATCH] scripts update
---
matlab/compile_results.m | 8 +-
matlab/compute/compute_fluxtube_growth_rate.m | 1 +
matlab/compute/mode_growth_meter.m | 2 +-
matlab/load/quick_gene_load.m | 4 +-
matlab/plot/create_film.m | 16 ++-
matlab/plot/photomaton.m | 13 +-
matlab/plot/plot_ballooning.m | 53 ++++++++
matlab/plot/poloidal_plot.m | 73 +++++++++++
matlab/plot/show_moments_spectrum.m | 25 +++-
matlab/process_field.m | 123 +++++++++---------
matlab/profiler.m | 59 +++++----
matlab/setup.m | 3 +-
matlab/write_fort90.m | 1 +
13 files changed, 272 insertions(+), 109 deletions(-)
create mode 100644 matlab/plot/plot_ballooning.m
create mode 100644 matlab/plot/poloidal_plot.m
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index b98ac6c1..8beff240 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -37,11 +37,12 @@ Ts5D_ = [];
Sipj_ = []; Sepj_ = [];
Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old;
Pi_max=0; Pe_max=0; Ji_max=0; Je_max=0;
+DATA.outfilenames = {};
+ii = 1;
while(CONTINUE)
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
- % Load results of simulation #JOBNUM
-% load_results
+ DATA.outfilenames{ii} = filename;
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Loading ',filename])
% Loading from output file
@@ -194,7 +195,8 @@ while(CONTINUE)
DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
-
+
+ ii = ii + 1;
JOBFOUND = JOBFOUND + 1;
LASTJOB = JOBNUM;
Pe_old = Pe_new; Je_old = Je_new;
diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index 60fb9437..b92ecbe9 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -47,6 +47,7 @@ if PLOT >0
xlim([min(linear_gr.ky) max(linear_gr.ky)]);
xlabel('$k_y$');
legend('show');
+ title(DATA.param_title);
if PLOT > 1
[YY,XX] = meshgrid(kys,t(its+1:ite));
figure
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 95af8e71..5dc00886 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -3,7 +3,7 @@ function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
t = OPTIONS.TIME;
-iz = 1;
+iz = OPTIONS.iz;
[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
FRAMES = zeros(size(OPTIONS.TIME));
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 3eddddf2..6d2dfa76 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -25,11 +25,11 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSGDK_0.0047_64x32.txt';
-fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
-% 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.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';
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index ed66da62..4c7ee671 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -8,8 +8,12 @@ else
T = DATA.Ts5D;
end
%% Processing
-toplot = process_field(DATA,OPTIONS);
-
+switch OPTIONS.PLAN
+ case 'RZ'
+ toplot = poloidal_plot(DATA,OPTIONS);
+ otherwise
+ toplot = process_field(DATA,OPTIONS);
+end
%%
FILENAME = [DATA.localdir,toplot.FILENAME,format];
XNAME = ['$',toplot.XNAME,'$'];
@@ -24,8 +28,8 @@ switch format
open(vidfile);
end
% Set colormap boundaries
-hmax = max(max(max(FIELD(:,:,FRAMES))));
-hmin = min(min(min(FIELD(:,:,FRAMES))));
+hmax = max(max(max(FIELD(:,:,:))));
+hmin = min(min(min(FIELD(:,:,:))));
scale = -1;
flag = 0;
if hmax == hmin
@@ -47,7 +51,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]);
end
in = 1;
nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
+ for n = 1:numel(FRAMES) % loop over selected frames
scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize
if ~strcmp(OPTIONS.PLAN,'sx')
if NORMALIZED
@@ -73,7 +77,7 @@ fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]);
else % show velocity distr.
contour(toplot.X,toplot.Y,FIELD(:,:,n)/scale,128);
end
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(FRAMES(n))))...
,', scaling = ',sprintf('%.1e',scale)]);
xlabel(XNAME); ylabel(YNAME); %colorbar;
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index 4a780f35..278fe024 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -1,7 +1,12 @@
function [ FIGURE ] = photomaton( DATA,OPTIONS )
%UNTITLED5 Summary of this function goes here
%% Processing
-toplot = process_field(DATA,OPTIONS);
+switch OPTIONS.PLAN
+ case 'RZ'
+ toplot = poloidal_plot(DATA,OPTIONS);
+ otherwise
+ toplot = process_field(DATA,OPTIONS);
+end
FNAME = toplot.FILENAME;
FRAMES = toplot.FRAMES;
Nframes= numel(FRAMES);
@@ -13,13 +18,13 @@ FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
for i_ = 1:numel(FRAMES)
subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
if ~strcmp(OPTIONS.PLAN,'kxky')
- scale = max(max(abs(toplot.FIELD(:,:,FRAMES(i_))))); % Scaling to normalize
+ scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize
else
scale = 1;
end
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');
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
if OPTIONS.AXISEQUAL
pbaspect(toplot.ASPECT)
end
@@ -31,7 +36,7 @@ FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
shading interp;
end
else
- contour(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale,128);
+ contour(toplot.X,toplot.Y,toplot.FIELD(:,:,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
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
new file mode 100644
index 00000000..f209aea5
--- /dev/null
+++ b/matlab/plot/plot_ballooning.m
@@ -0,0 +1,53 @@
+function [FIG] = plot_ballooning(data,options)
+
+ [~,it] = min(abs(data.Ts3D - options.time_2_plot));
+ phi_real=(real(data.PHI(:,:,:,it)));
+ phi_imag=(imag(data.PHI(:,:,:,it)));
+ % Apply baollooning tranform
+ for iky=2
+ dims = size(phi_real);
+ phib_real = zeros( dims(1)*dims(3) ,1);
+ phib_imag= phib_real;
+ b_angle = phib_real;
+
+ midpoint = floor((dims(1)*dims(3) )/2)+1;
+ for ip =1: dims(1)
+ start_ = (ip -1)*dims(3) +1;
+ end_ = ip*dims(3);
+ phib_real(start_:end_) = phi_real(ip,iky,:);
+ phib_imag(start_:end_) = phi_imag(ip,iky, :);
+ end
+
+ % Define ballooning angle
+ Nkx = numel(data.kx)-1; coordz = data.z;
+ idx = -Nkx:1:Nkx;
+ for ip =1: dims(1)
+ for iz=1:dims(3)
+ ii = dims(3)*(ip -1) + iz;
+ b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
+ end
+ end
+
+ phib = phib_real(:) + 1i * phib_imag(:);
+ % normalize real and imaginary parts at chi =0
+ if options.normalized
+ [~,idxLFS] = min(abs(b_angle -0));
+ normalization = phib( idxLFS);
+ else
+ normalization = 1;
+ end
+ phib_norm = phib / normalization ;
+ phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
+ phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
+
+
+ FIG.fig = figure; hold all;
+ plot(b_angle/pi, phib_real_norm,'-b');
+ plot(b_angle/pi, phib_imag_norm ,'-r');
+ plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k');
+ legend('real','imag','norm')
+ xlabel('$\chi / \pi$')
+ ylabel('$\phi_B (\chi)$');
+ title(['HeLaZ ballooning, t=',num2str(data.Ts3D(it))]);
+ end
+end
diff --git a/matlab/plot/poloidal_plot.m b/matlab/plot/poloidal_plot.m
new file mode 100644
index 00000000..66775a91
--- /dev/null
+++ b/matlab/plot/poloidal_plot.m
@@ -0,0 +1,73 @@
+function [TOPLOT] = poloidal_plot(DATA,OPTIONS)
+sq = @(x) squeeze(x);
+%% Geometry
+r0 = DATA.a;
+q0 = DATA.Q0;
+Cy = r0/q0;
+Ly = max(DATA.y) - min(DATA.y);
+n0 = Cy*2*pi/Ly;
+z_ = DATA.z;
+%% grid sizes
+nkx = DATA.Nkx; nky = DATA.Nky; nz = DATA.Nz;
+nx = DATA.Nx; ny = DATA.Ny;
+%% Time and frames
+FRAMES = zeros(size(OPTIONS.TIME));
+for i = 1:numel(OPTIONS.TIME)
+ [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
+end
+FRAMES = unique(FRAMES);
+
+%% Field to plot
+switch OPTIONS.NAME
+ case '\phi'
+ FIELD = DATA.PHI;
+ NAME = 'phi';
+ case 'n_i'
+ FIELD = DATA.DENS_I;
+ NAME = 'ni';
+end
+
+%% Build poloidal plot (tor angle = 0)
+FRZ = zeros(DATA.Nx,DATA.Nz, numel(FRAMES));
+w_ = zeros(nky,nz);
+for iz = 1:nz
+ for iky = 1:ny
+ w_(iky,iz) = exp(1i*iky*n0*(q0*z_(iz)));
+ end
+end
+
+
+for it = 1:numel(FRAMES)
+ field_ = squeeze(FIELD(:,:,:,FRAMES(it)));
+ spectrumKxKyZ = zeros(nkx,nky,nz);
+ spectrumKxKyZ(1:nkx,:,:) = field_;
+ spectrumKxKyZ((nkx+1):nx,1,:) = conj(field_(nkx:-1:2,1,:));
+ spectrumKxKyZ((nkx+1):nx,2:ny,:) = conj(field_(nkx:-1:2,ny:-1:2,:));
+
+ Fxkyz = ny*ifft(spectrumKxKyZ,[],1);
+ for ix = 1:DATA.Nx
+ FRZ(ix,:,it) = sq(sum(sq(Fxkyz(ix,:,:)).*w_(:,:),1));
+ end
+end
+FRZ = real(FRZ);
+%grid
+[Y,X] = meshgrid(DATA.z,DATA.x);
+X__ = (X+DATA.a).*cos(Y);
+Y__ = (X+DATA.a).*sin(Y);
+X = X__;
+Y = Y__;
+%% output
+
+TOPLOT.FIELD = FRZ;
+TOPLOT.FRAMES = FRAMES;
+TOPLOT.X = (X);
+TOPLOT.Y = (Y);
+TOPLOT.FIELDNAME = NAME;
+TOPLOT.XNAME = '$R$';
+TOPLOT.YNAME = '$Z$';
+TOPLOT.FILENAME = [NAME,'_poloidal_phieq0'];
+TOPLOT.DIMENSIONS= [100, 100, 500, 500];
+TOPLOT.ASPECT = [1 1 1];
+TOPLOT.INTERP = OPTIONS.INTERP;
+
+end
\ No newline at end of file
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 5dc78fae..6e80a71f 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -6,10 +6,12 @@ Ji = DATA.Ji;
Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
Nipj = squeeze(Nipj);
+if DATA.K_E
Pe = DATA.Pe;
Je = DATA.Je;
Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
Nepj = squeeze(Nepj);
+end
FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS];
set(gcf, 'Position', [100 50 1000 400])
@@ -20,22 +22,27 @@ if ~OPTIONS.ST
[~,it0] = min(abs(t0-DATA.Ts5D));
[~,it1] = min(abs(t1-DATA.Ts5D));
Nipj = mean(Nipj(:,:,it0:it1),3);
+ if DATA.K_E
Nepj = mean(Nepj(:,:,it0:it1),3);
+ end
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);
+
+ if DATA.K_E
+ Nepj = squeeze(Nepj);
ymine = min(Nepj); ymaxe = max(Nepj);
- ymin = min([ymini ymine]);
ymax = max([ymaxi ymaxe]);
-
-
+ ymin = min([ymini ymine]);
+ end
+ if DATA.K_E
subplot(121)
+ end
if ~OPTIONS.P2J
for ij = 1:numel(Ji)
name = ['$j=',num2str(Ji(ij)),'$'];
@@ -52,7 +59,8 @@ if ~OPTIONS.ST
ylabel(['$\sum_{kx,ky}|N_i^{pj}|$']);
legend('show');
title([TITLE,' He-La ion spectrum']);
-
+
+ if DATA.K_E
subplot(122)
if ~OPTIONS.P2J
for ij = 1:numel(Je)
@@ -70,6 +78,7 @@ if ~OPTIONS.ST
ylabel(['$\sum_{kx,ky}|N_e^{pj}|$']);
legend('show');
title([TITLE,' He-La elec. spectrum']);
+ end
else
[JJ,PP] = meshgrid(Ji,Pi);
P2Ji = PP + 2*JJ;
@@ -91,6 +100,7 @@ else
for ip2j = 1:numel(p2ji)
Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
end
+ if DATA.K_E
% same for electrons!!
[JJ,PP] = meshgrid(Je,Pe);
P2Je = PP + 2*JJ;
@@ -112,6 +122,7 @@ else
for ip2j = 1:numel(p2ji)
Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
end
+ end
% plots
% scaling
% plt = @(x,ip2j) x./max(max(x));
@@ -120,13 +131,16 @@ else
else
plt = @(x,ip2j) x;
end
+ if DATA.K_E
subplot(2,1,1)
+ end
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}$')
+ if DATA.K_E
subplot(2,1,2)
imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ji)));
set(gca,'YDir','normal')
@@ -135,6 +149,7 @@ else
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/process_field.m b/matlab/process_field.m
index f19f1b2d..5762fc84 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -7,10 +7,11 @@ FRAMES = zeros(size(OPTIONS.TIME));
for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
end
+
%% Setup the plot geometry
[KY,~] = meshgrid(DATA.ky,DATA.kx);
directions = {'x','y','z'};
-Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(OPTIONS.TIME);
+Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
@@ -134,8 +135,8 @@ switch OPTIONS.NAME
case '\phi'
NAME = 'phi';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.PHI(:,:,:,it)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -144,9 +145,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -154,8 +155,8 @@ switch OPTIONS.NAME
case 'N_i^{00}'
NAME = 'Ni00';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.Ni00(:,:,:,it)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.Ni00(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -164,9 +165,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -174,8 +175,8 @@ switch OPTIONS.NAME
case 'n_e'
NAME = 'ne';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.DENS_E(:,:,:,it)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -184,9 +185,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -195,9 +196,9 @@ switch OPTIONS.NAME
NAME = 'k2ne';
[KY, KX] = meshgrid(DATA.ky, DATA.kx);
if COMPDIM == 3
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it)));
+ tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(process(tmp));
end
@@ -207,9 +208,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -217,8 +218,8 @@ switch OPTIONS.NAME
case 'n_e^{NZ}'
NAME = 'neNZ';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.DENS_E(:,:,:,it).*(KY~=0)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -227,9 +228,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it).*(KY~=0)));
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -237,8 +238,8 @@ switch OPTIONS.NAME
case 'n_i'
NAME = 'ni';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.DENS_I(:,:,:,it)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -247,9 +248,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -257,8 +258,8 @@ switch OPTIONS.NAME
case 'T_i'
NAME = 'Ti';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.TEMP_I(:,:,:,it)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.TEMP_I(:,:,:,FRAMES(it))));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -267,9 +268,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,it)));
+ tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,FRAMES(it))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -277,8 +278,8 @@ switch OPTIONS.NAME
case 'n_i^{NZ}'
NAME = 'niNZ';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.DENS_I(:,:,:,it).*(KY~=0)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -287,9 +288,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it).*(KY~=0)));
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -297,8 +298,8 @@ switch OPTIONS.NAME
case '\phi^{NZ}'
NAME = 'phiNZ';
if COMPDIM == 3
- for it = FRAMES
- tmp = squeeze(compr(DATA.PHI(:,:,:,it).*(KY~=0)));
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it)).*(KY~=0)));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
@@ -307,9 +308,9 @@ switch OPTIONS.NAME
else
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
end
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it).*(KY~=0)));
+ tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it)).*(KY~=0)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -318,21 +319,21 @@ switch OPTIONS.NAME
NAME = 'vy';
[~, 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 it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
end
end
if REALP % plot in real space
- for it = FRAMES
+ 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 = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,it))));
+ tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -341,21 +342,21 @@ switch OPTIONS.NAME
NAME = 'vx';
[KY, ~] = 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 it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
end
end
if REALP % plot in real space
- for it = FRAMES
+ 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 = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,it))));
+ tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -364,21 +365,21 @@ switch OPTIONS.NAME
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 it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny));
end
end
if REALP % plot in real space
- for it = FRAMES
+ 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 = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,it))));
+ tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
@@ -387,14 +388,14 @@ switch OPTIONS.NAME
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 it = 1:numel(FRAMES)
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)));
+ 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)));
end
end
if REALP % plot in real space
- for it = FRAMES
+ for it = 1:numel(FRAMES)
FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
end
else % Plot spectrum
@@ -402,7 +403,7 @@ switch OPTIONS.NAME
shift_x = @(x) fftshift(x,2);
shift_y = @(x) fftshift(x,2);
tmp = zeros(DATA.Nx,DATA.Ny,Nz);
- for it = FRAMES
+ for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
end
@@ -417,8 +418,8 @@ switch OPTIONS.NAME
dit_ = 1;%ceil((it1_-it0_+1)/10);
FRAMES = it0_:dit_:it1_;
iz = 1;
- for it = FRAMES
- OPTIONS.T = DATA.Ts5D(it);
+ for it = 1:numel(FRAMES)
+ OPTIONS.T = DATA.Ts5D(FRAMES(it));
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
case 'f_e'
@@ -429,8 +430,8 @@ switch OPTIONS.NAME
dit_ = 1;%ceil((it1_-it0_+1)/10);
FRAMES = it0_:dit_:it1_;
iz = 1;
- for it = FRAMES
- OPTIONS.T = DATA.Ts5D(it);
+ for it = 1:numel(FRAMES)
+ OPTIONS.T = DATA.Ts5D(FRAMES(it));
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
end
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 2988e22c..a8184e13 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -1,31 +1,34 @@
%% load profiling
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
-% filename_ = ['/misc/HeLaZ_outputs',filename(3:end)];
-filename_ = ['/home/ahoffman/HeLaZ',filename(3:end)];
-CPUTIME = double(h5readatt(filename_,'/data/input','cpu_time'));
-DT_SIM = h5readatt(filename_,'/data/input','dt');
+% outfilename = ['/misc/HeLaZ_outputs',filename(3:end)];
+outfilename = data.outfilenames{end};
+CPUTIME = double(h5readatt(outfilename,'/data/input','cpu_time'));
+DT_SIM = h5readatt(outfilename,'/data/input','dt');
-rhs_Tc = h5read(filename_,'/profiler/Tc_rhs');
-adv_field_Tc = h5read(filename_,'/profiler/Tc_adv_field');
-ghost_Tc = h5read(filename_,'/profiler/Tc_ghost');
-clos_Tc = h5read(filename_,'/profiler/Tc_clos');
-coll_Tc = h5read(filename_,'/profiler/Tc_coll');
-poisson_Tc = h5read(filename_,'/profiler/Tc_poisson');
-Sapj_Tc = h5read(filename_,'/profiler/Tc_Sapj');
-checkfield_Tc= h5read(filename_,'/profiler/Tc_checkfield');
-diag_Tc = h5read(filename_,'/profiler/Tc_diag');
-step_Tc = h5read(filename_,'/profiler/Tc_step');
-Ts0D = h5read(filename_,'/profiler/time');
+rhs_Tc = h5read(outfilename,'/profiler/Tc_rhs');
+adv_field_Tc = h5read(outfilename,'/profiler/Tc_adv_field');
+ghost_Tc = h5read(outfilename,'/profiler/Tc_ghost');
+clos_Tc = h5read(outfilename,'/profiler/Tc_clos');
+coll_Tc = h5read(outfilename,'/profiler/Tc_coll');
+poisson_Tc = h5read(outfilename,'/profiler/Tc_poisson');
+Sapj_Tc = h5read(outfilename,'/profiler/Tc_Sapj');
+checkfield_Tc= h5read(outfilename,'/profiler/Tc_checkfield');
+diag_Tc = h5read(outfilename,'/profiler/Tc_diag');
+process_Tc = h5read(outfilename,'/profiler/Tc_process');
+step_Tc = h5read(outfilename,'/profiler/Tc_step');
+Ts0D = h5read(outfilename,'/profiler/time');
+
+N_T = 11;
missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ...
- -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc;
+ -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc-process_Tc;
total_Tc = step_Tc;
TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(ghost_Tc);...
diff(clos_Tc); diff(coll_Tc); diff(poisson_Tc); diff(Sapj_Tc); ...
- diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
-TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/10,10]);
+ diff(checkfield_Tc); diff(diag_Tc); diff(process_Tc); diff(missing_Tc)];
+TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]);
TIME_PER_STEP = sum(TIME_PER_FCT,2);
TIME_PER_CPU = trapz(Ts0D(2:end),TIME_PER_STEP);
@@ -38,6 +41,7 @@ coll_Ta = mean(diff(coll_Tc));
poisson_Ta = mean(diff(poisson_Tc));
Sapj_Ta = mean(diff(Sapj_Tc));
checkfield_Ta = mean(diff(checkfield_Tc));
+process_Ta = mean(diff(process_Tc));
diag_Ta = mean(diff(diag_Tc));
NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM;
@@ -46,23 +50,24 @@ NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM;
if 1
%% Area plot
fig = figure;
-
+% colors = rand(N_T,3);
+colors = lines(N_T);
p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none');
-for i = 1:10; p1(i).FaceColor = rand(1,3); end;
-legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
+for i = 1:N_T; p1(i).FaceColor = colors(i,:); end;
+legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
xlim([Ts0D(2),Ts0D(end)]);
title(sprintf('Proc. 1, total sim. time ~%.0f [h]',CPUTIME/3600))
hold on
FIGNAME = 'profiler';
-save_figure
+% save_figure
else
%% Normalized Area plot
fig = figure;
p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none', 'FaceColor','flat');
-for i = 1:10; p1(i).FaceColor = rand(1,3); end;
+for i = 1:N_T; p1(i).FaceColor = rand(1,3); end;
legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
@@ -73,7 +78,7 @@ ylabel('Step Comp. Time [s]')
ylim([0,1.1*max(diff(total_Tc))])
set(gca,'ycolor','r')
FIGNAME = 'profiler';
-save_figure
+% save_figure
end
if 0
@@ -86,11 +91,13 @@ histogram(diff(clos_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(coll_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(poisson_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(Sapj_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(process_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(checkfield_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(diag_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
grid on;
-legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
+legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Process','Check+sym', 'Diagnos.', 'Missing')
xlabel('Step Comp. Time [s]'); ylabel('')
+set(gca,'Xscale','log')
FIGNAME = 'profiler';
-save_figure
+% save_figure
end
\ No newline at end of file
diff --git a/matlab/setup.m b/matlab/setup.m
index 00269b32..5233fad4 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -1,4 +1,4 @@
-%% ________________________________________________________________________
+%% _______________________________________________________________________
SIMDIR = ['../results/',SIMID,'/'];
% Grid parameters
GRID.pmaxe = PMAXE; % Electron Hermite moments
@@ -10,6 +10,7 @@ GRID.Lx = LX; % x length
GRID.Ny = NY; % y ''
GRID.Ly = LY; % y ''
GRID.Nz = NZ; % z resolution
+GRID.Npol = NPOL; % z resolution
if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
% Geometry
GEOM.geom = ['''',GEOMETRY,''''];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index d35641f0..26439459 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -20,6 +20,7 @@ fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']);
fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']);
fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
+fprintf(fid,[' Npol = ', num2str(GRID.Npol),'\n']);
fprintf(fid,[' SG = ', GRID.SG,'\n']);
fprintf(fid,'/\n');
--
GitLab