From 843bee0072643473c3765c123f8fc055949d6a9a Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 10 Oct 2022 11:12:12 +0200
Subject: [PATCH] fix some matlab scripts
---
matlab/evolve_tracers.m | 18 +++++------
matlab/extract_fig_data.m | 27 +++++++++++-----
matlab/load/load_grid_data.m | 5 ++-
.../plot_radial_transport_and_spacetime.m | 8 +++--
matlab/process_field.m | 31 ++++++++++++++++---
5 files changed, 64 insertions(+), 25 deletions(-)
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 952a87ad..a6a738ce 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,19 +1,19 @@
% Options
-SHOW_FILM = 0;
-field2plot ='Ni00';
+SHOW_FILM = 1;
+field2plot ='phi';
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME = 200; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME = 40; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin = 200;
-dt_ = 0.1;
+Tfin = 100;
+dt_ = 0.01;
Nstep = ceil(Tfin/dt_);
% Init tracers
Np = 20; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
-
-Na = 1/dt_; %length of trace
+dimmed = 0; % to dimm the colormap in the background (infty = white, 0 normal color)
+Na = 10/dt_; %length of trace
Traj_x = zeros(Np,Nstep);
Traj_y = zeros(Np,Nstep);
@@ -224,8 +224,8 @@ while(t_<Tfin && it <= Nstep)
end
scale = max(max(abs(F2P))); % Scaling to normalize
pclr = pcolor(XX_,YY_,F2P/scale);
- colormap(bluewhitered); caxis([-1 1]);
- set(pclr, 'edgecolor','none'); hold on; caxis([-2,2]); shading interp
+ caxis([-1 1]); colormap(bluewhitered); colorbar;
+ set(pclr, 'edgecolor','none'); hold on; caxis([-1+dimmed,1+dimmed]); shading interp
for ip = 1:Np
ia0 = max(1,it-Na);
plot(Traj_x(ip,ia0:it),Traj_y(ip,ia0:it),'.','Color',color(ip,:)); hold on
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index e7a98b9d..7cd704c8 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -1,3 +1,11 @@
+% tw = [500 1000];
+% tw = [1000 1500];
+% tw = [2000 3000];
+% tw = [3000 4000];
+% tw = [4000 4500];
+% tw = [4500 5000];
+tw = [0 6500];
+
fig = gcf;
axObjs = fig.Children;
dataObjs = axObjs.Children;
@@ -7,18 +15,21 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData];
end
-n0 = 1;
-n1 = numel(X_);
-figure;
-mvm = @(x) movmean(x,1);
+% n0 = 1;
+% n1 = numel(X_);
+[~,n0] = min(abs(X_-tw(1)));
+[~,n1] = min(abs(X_-tw(2)));
+mvm = @(x) movmean(x,50);
shift = X_(n0);
% shift = 0;
-% plot(X_(n0:end),Y_(n0:end));
+
+figure;
+plot(X_(n0:end),Y_(n0:end));
plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
-t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
-avg= mean(Y_(t0:t1)); dev = std(Y_(t0:t1));
- disp(['AVG =',sprintf('%2.2f',avg),'+-',sprintf('%2.2f',dev)]);
+% t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
+avg= mean(Y_(n0:n1)); dev = std(Y_(n0:n1));
+ disp(['AVG =',sprintf('%4.4f',avg),'+-',sprintf('%4.4f',dev)]);
%
% n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
% avg_ = mean(Y_(n1:n2));
diff --git a/matlab/load/load_grid_data.m b/matlab/load/load_grid_data.m
index be964241..940381ec 100644
--- a/matlab/load/load_grid_data.m
+++ b/matlab/load/load_grid_data.m
@@ -1,4 +1,4 @@
-function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data_3D( filename )
+function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data( filename )
%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ
pe = h5read(filename,'/data/grid/coordp_e');
je = h5read(filename,'/data/grid/coordj_e');
@@ -7,4 +7,7 @@ function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data_3D( filename )
kx = h5read(filename,'/data/grid/coordkx');
ky = h5read(filename,'/data/grid/coordky');
z = h5read(filename,'/data/grid/coordz');
+ if (numel(z) == 1)
+ z = 0;
+ end
end
\ No newline at end of file
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index de18d2d5..5ea5cbc5 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -10,8 +10,8 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
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;
- disp(['G_x=',sprintf('%2.2f',Gx_infty_avg),'+-',sprintf('%2.2f',Gx_infty_std)]);
- disp(['Q_x=',sprintf('%2.2f',Qx_infty_avg),'+-',sprintf('%2.2f',Qx_infty_std)]);
+ disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
+% disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
ikzf = min([ikzf,DATA.Nky]);
@@ -27,7 +27,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
end
Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee);
- disp(['Q_avg=',sprintf('%2.2f',Qx_avg),'+-',sprintf('%2.2f',Qx_err)]);
+% disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
%% computations
% Compute Gamma from ifft matlab
@@ -78,7 +78,9 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
ylabel('$Q_x$')
+ if(~isnan(Qx_infty_avg))
ylim([0,5*abs(Qx_infty_avg)]);
+ end
xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
diff --git a/matlab/process_field.m b/matlab/process_field.m
index ed7fe9e3..0804e892 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -52,8 +52,8 @@ 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];
+DIMENSIONS = [100, 100, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
% Polar grid
POLARNAME = [];
if POLARPLOT
@@ -64,7 +64,7 @@ if POLARPLOT
Y = Y__;
XNAME='X';
YNAME='Z';
- DIMENSIONS = [100, 100, 500, 500];
+ DIMENSIONS = [100, 100, OPTIONS.RESOLUTION, OPTIONS.RESOLUTION];
ASPECT = [1,1,1];
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
@@ -407,8 +407,8 @@ switch OPTIONS.NAME
end
case 's_y'
NAME = 'sy';
- [~, 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(ifourier_GENE(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
@@ -428,6 +428,29 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case '\omega_z'
+ NAME = 'vorticity';
+ [KX, KY] = meshgrid(DATA.kx, DATA.ky);
+ wE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
+ for it = 1:numel(FRAMES)
+ for iz = 1:DATA.Nz
+ wE(:,:,iz,it) = real(ifourier_GENE(-(KX.^2+KY.^2).*(DATA.PHI(:,:,iz,FRAMES(it)))));
+ end
+ end
+ if REALP % plot in real space
+ for it = 1:numel(FRAMES)
+ FIELD(:,:,it) = squeeze(compr(wE(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
+ for it = 1:FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*(DATA.PHI(:,:,iz,FRAMES(it)))));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
case '\Gamma_x'
NAME = 'Gx';
[~, KY] = meshgrid(DATA.kx, DATA.ky);
--
GitLab