From f4ce4a1ac1639119a51da309f7cc90b66616ea61 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 8 Nov 2021 11:37:26 +0100
Subject: [PATCH] update scripts
---
matlab/compile_results.m | 143 +++++++----
matlab/{create_gif.m => create_film.m} | 65 +++--
matlab/create_mov.m | 69 ------
matlab/load_marconi.m | 3 +-
matlab/load_params.m | 103 ++++----
matlab/photomaton.m | 24 ++
.../plots/plot_radial_transport_and_shear.m | 104 ++++----
matlab/post_processing.m | 20 --
matlab/process_field.m | 196 +++++++++++++++
matlab/save_figure.m | 12 +-
matlab/show_geometry.m | 96 ++++++++
wk/analysis_3D.m | 230 ++++--------------
wk/linear_1D_entropy_mode.m | 13 +-
wk/local_run.m | 14 +-
wk/marconi_run.m | 46 ++--
wk/photomaton.m | 59 -----
wk/shearless_linear_fluxtube.m | 77 ++++--
17 files changed, 714 insertions(+), 560 deletions(-)
rename matlab/{create_gif.m => create_film.m} (52%)
delete mode 100644 matlab/create_mov.m
create mode 100644 matlab/photomaton.m
create mode 100644 matlab/process_field.m
create mode 100644 matlab/show_geometry.m
delete mode 100644 wk/photomaton.m
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index b4c3baeb..dc515acd 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -1,12 +1,14 @@
+function [DATA] = compile_results(DIRECTORY,JOBNUMMIN,JOBNUMMAX)
+DATA = {};
CONTINUE = 1;
JOBNUM = JOBNUMMIN; JOBFOUND = 0;
-TJOB_SE = []; % Start and end times of jobs
-NU_EVOL = []; % evolution of parameter nu between jobs
-CO_EVOL = []; % evolution of CO
-MU_EVOL = []; % evolution of parameter mu between jobs
-K_N_EVOL = []; %
-L_EVOL = []; %
-DT_EVOL = []; %
+DATA.TJOB_SE = []; % Start and end times of jobs
+DATA.NU_EVOL = []; % evolution of parameter nu between jobs
+DATA.CO_EVOL = []; % evolution of CO
+DATA.MU_EVOL = []; % evolution of parameter mu between jobs
+DATA.K_N_EVOL = []; %
+DATA.L_EVOL = []; %
+DATA.DT_EVOL = []; %
% FIELDS
Nipj_ = []; Nepj_ = [];
Ni00_ = []; Ne00_ = [];
@@ -25,10 +27,27 @@ 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;
while(CONTINUE)
- filename = sprintf([BASIC.MISCDIR,'outputs_%.2d.h5'],JOBNUM);
+ filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
% Load results of simulation #JOBNUM
- load_results
+% load_results
+ %% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ disp(['Loading ',filename])
+ % Loading from output file
+ CPUTIME = h5readatt(filename,'/data/input','cpu_time');
+ DT_SIM = h5readatt(filename,'/data/input','dt');
+ [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
+ W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
+ W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
+ W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
+ W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
+ W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
+ W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
+ W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
+ W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
+% KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
+ KIN_E = 1;
+
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
Pi_new= numel(Pi); Ji_new= numel(Ji);
@@ -74,50 +93,54 @@ while(CONTINUE)
end
end
- if W_GAMMA || W_HF
- Ts0D_ = cat(1,Ts0D_,Ts0D);
- end
if W_GAMMA
- GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI);
- PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI);
+ [ GGAMMA_RI, Ts0D, ~] = load_0D_data(filename, 'gflux_ri');
+ PGAMMA_RI = load_0D_data(filename, 'pflux_ri');
+ GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); clear GGAMMA_RI
+ PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); clear PGAMMA_RI
end
if W_HF
- HFLUX_ = cat(1,HFLUX_,HFLUX_X);
+ [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x');
+ HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X
end
- if W_PHI || W_NA00
- Ts3D_ = cat(1,Ts3D_,Ts3D);
- end
if W_PHI
- PHI_ = cat(4,PHI_,PHI);
+ [ PHI, Ts3D, ~] = load_3D_data(filename, 'phi');
+ PHI_ = cat(4,PHI_,PHI); clear PHI
end
if W_NA00
if KIN_E
- Ne00_ = cat(4,Ne00_,Ne00);
+ Ne00 = load_3D_data(filename, 'Ne00');
+ Ne00_ = cat(4,Ne00_,Ne00); clear Ne00
end
- Ni00_ = cat(4,Ni00_,Ni00);
+ [Ni00, Ts3D, ~] = load_3D_data(filename, 'Ni00');
+ Ni00_ = cat(4,Ni00_,Ni00); clear Ni00
end
if W_DENS
if KIN_E
- DENS_E_ = cat(4,DENS_E_,DENS_E);
+ [DENS_E, Ts3D, ~] = load_3D_data(filename, 'dens_e');
+ DENS_E_ = cat(4,DENS_E_,DENS_E); clear DENS_E
end
- DENS_I_ = cat(4,DENS_I_,DENS_I);
+ [DENS_I, Ts3D, ~] = load_3D_data(filename, 'dens_i');
+ DENS_I_ = cat(4,DENS_I_,DENS_I); clear DENS_I
end
if W_TEMP
if KIN_E
- TEMP_E_ = cat(4,TEMP_E_,TEMP_E);
+ [TEMP_E, Ts3D, ~] = load_3D_data(filename, 'temp_e');
+ TEMP_E_ = cat(4,TEMP_E_,TEMP_E); clear TEMP_E
end
- TEMP_I_ = cat(4,TEMP_I_,TEMP_I);
- end
- if W_NAPJ || W_SAPJ
- Ts5D_ = cat(1,Ts5D_,Ts5D);
+ [TEMP_I, Ts3D, ~] = load_3D_data(filename, 'temp_i');
+ TEMP_I_ = cat(4,TEMP_I_,TEMP_I); clear TEMP_I
end
+
if W_NAPJ
- Nipj_ = cat(6,Nipj_,Nipj);
+ [Nipj, Ts5D, ~] = load_5D_data(filename, 'moments_i');
+ Nipj_ = cat(6,Nipj_,Nipj); clear Nipj
if KIN_E
- Nepj_ = cat(6,Nepj_,Nepj);
+ Nepj = load_5D_data(filename, 'moments_e');
+ Nepj_ = cat(6,Nepj_,Nepj); clear Nepj
end
end
if W_SAPJ
@@ -126,16 +149,19 @@ while(CONTINUE)
Sepj_ = cat(6,Sepj_,Sepj);
end
end
+ Ts0D_ = cat(1,Ts0D_,Ts0D);
+ Ts3D_ = cat(1,Ts3D_,Ts3D);
+ Ts5D_ = cat(1,Ts5D_,Ts5D);
% Evolution of simulation parameters
load_params
- TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
- NU_EVOL = [NU_EVOL NU NU];
- CO_EVOL = [CO_EVOL CO CO];
- MU_EVOL = [MU_EVOL MU MU];
- K_N_EVOL = [K_N_EVOL K_N K_N];
- L_EVOL = [L_EVOL L L];
- DT_EVOL = [DT_EVOL DT_SIM DT_SIM];
+ DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)];
+ DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU];
+ DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
+ DATA.MU_EVOL = [DATA.MU_EVOL DATA.MU DATA.MU];
+ 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];
JOBFOUND = JOBFOUND + 1;
LASTJOB = JOBNUM;
@@ -147,13 +173,40 @@ while(CONTINUE)
end
JOBNUM = JOBNUM + 1;
end
-GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; HFLUX_X = HFLUX_; Ts0D = Ts0D_;
-Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
-Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts3D = Ts3D_;
-DENS_E = DENS_E_; DENS_I = DENS_I_; TEMP_E = TEMP_E_; TEMP_I = TEMP_I_;
-clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_ HFLUX_
-Sipj = Sipj_; Sepj = Sepj_;
-clear Sipj_ Sepj_
+%% Build grids
+Nkx = numel(kx); Nky = numel(ky);
+[KY,KX] = meshgrid(ky,kx);
+Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
+dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
+KPERP2 = KY.^2+KX.^2;
+[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
+
+Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z);
+Lx = 2*pi/dkx; Ly = 2*pi/dky;
+dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz;
+x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
+y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
+%% Add everything in output structure
+% Fields
+DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_;
+DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
+if(KIN_E)
+DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_;
+end
+DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_;
+DATA.PHI = PHI_;
+% grids
+DATA.Pe = Pe; DATA.Pi = Pi;
+DATA.kx = kx; DATA.ky = ky; DATA.z = z;
+DATA.x = x; DATA.y = y;
+DATA.ikx0 = ikx0; DATA.iky0 = iky0;
+DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
+DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji);
+DATA.dir = DIRECTORY;
+DATA.localdir = ['..',DIRECTORY(20:end)];
+
JOBNUM = LASTJOB;
-filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
+
+filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+end
\ No newline at end of file
diff --git a/matlab/create_gif.m b/matlab/create_film.m
similarity index 52%
rename from matlab/create_gif.m
rename to matlab/create_film.m
index 1d0d5717..b48cf206 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_film.m
@@ -1,11 +1,24 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-XNAME = latexize(XNAME);
-YNAME = latexize(YNAME);
-FIELDNAME = latexize(FIELDNAME);
+function create_film(DATA,OPTIONS,format)
+%% Plot options
+FPS = 30; DELAY = 1/FPS;
+INTERP = OPTIONS.INTERP; BWR = 1; NORMALIZED = 1;
+T = DATA.Ts3D;
+%% Processing
+toplot = process_field(DATA,OPTIONS);
+%%
+FILENAME = [toplot.FILENAME,format];
+XNAME = ['$',toplot.XNAME,'$'];
+YNAME = ['$',toplot.YNAME,'$'];
+FIELDNAME = ['$',toplot.FIELDNAME,'$'];
+FIELD = toplot.FIELD; X = toplot.X; Y = toplot.Y;
+FRAMES = toplot.FRAMES;
+switch format
+ case '.avi'
+ vidfile = VideoWriter(FILENAME,'Uncompressed AVI');
+ vidfile.FrameRate = FPS;
+ open(vidfile);
+end
% Set colormap boundaries
hmax = max(max(max(FIELD(:,:,FRAMES))));
hmin = min(min(min(FIELD(:,:,FRAMES))));
@@ -15,7 +28,7 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-fig = figure('Color','white','Position', [100, 100, 400, 400]);
+fig = figure('Color','white','Position', toplot.DIMENSIONS);
pcolor(X,Y,FIELD(:,:,1)); % to set up
if BWR
colormap(bluewhitered)
@@ -46,7 +59,7 @@ fig = figure('Color','white','Position', [100, 100, 400, 400]);
if INTERP
shading interp;
end
- set(pclr, 'edgecolor','none'); axis square;
+ set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT);
if BWR
colormap(bluewhitered)
end
@@ -54,14 +67,22 @@ fig = figure('Color','white','Position', [100, 100, 400, 400]);
drawnow
% Capture the plot as an image
frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,32);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
+ switch format
+ case '.gif'
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,32);
+ % Write to the GIF File
+ if in == 1
+ imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ case '.avi'
+ writeVideo(vidfile,frame);
+ otherwise
+ disp('Unknown format');
+ break
+ end
% terminal info
while nbytes > 0
fprintf('\b')
@@ -71,6 +92,12 @@ fig = figure('Color','white','Position', [100, 100, 400, 400]);
in = in + 1;
end
disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
+ switch format
+ case '.gif'
+ disp(['Gif saved @ : ',FILENAME])
+ case '.avi'
+ disp(['Video saved @ : ',FILENAME])
+ close(vidfile);
+ end
+end
end
-
diff --git a/matlab/create_mov.m b/matlab/create_mov.m
deleted file mode 100644
index 5f84d65d..00000000
--- a/matlab/create_mov.m
+++ /dev/null
@@ -1,69 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME];
-XNAME = latexize(XNAME);
-YNAME = latexize(YNAME);
-FIELDNAME = latexize(FIELDNAME);
-
-vidfile = VideoWriter(GIFNAME,'Uncompressed AVI');
-vidfile.FrameRate = FPS;
-open(vidfile);
-% Set colormap boundaries
-hmax = max(max(max(FIELD(:,:,FRAMES))));
-hmin = min(min(min(FIELD(:,:,FRAMES))));
-scale = -1;
-flag = 0;
-if hmax == hmin
- disp('Warning : h = hmin = hmax = const')
-else
-% Setup figure frame
-figure('Color','white','Position', [100, 100, 500, 500]);
- pcolor(X,Y,FIELD(:,:,1)); % to set up
-% colormap gray
- colormap(bluewhitered)
- axis tight manual % this ensures that getframe() returns a consistent size
- if INTERP
- shading interp;
- end
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
- scale = max(max(abs(FIELD(:,:,n))));
- pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot
- if INTERP
- shading interp;
- end
- if NORMALIZED
- pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
- caxis([-1,1]);
- else
- pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
- if CONST_CMAP
- caxis([-1,1]*max(abs([hmin hmax]))); % const color map
- else
- caxis([-1,1]*scale); % adaptive color map
- end
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scale = ',sprintf('%.1e',scale)]);
- end
- set(pclr, 'edgecolor','none'); axis square;
- xlabel(XNAME); ylabel(YNAME); %colorbar;
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale)]);
- drawnow
- % Capture the plot as an image
- frame = getframe(gcf);
- writeVideo(vidfile,frame);
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
- in = in + 1;
- end
- disp(' ')
- disp(['Video saved @ : ',GIFNAME])
- close(vidfile);
-end
-
diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m
index faee4369..2fb78159 100644
--- a/matlab/load_marconi.m
+++ b/matlab/load_marconi.m
@@ -16,7 +16,8 @@ function [ RESDIR ] = load_marconi( outfilename )
% disp(CMD);
% system(CMD);
% SCP the output file from marconi to misc folder of SPCPC
- CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
+% CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
+ CMD = ['rsync ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
disp(CMD);
system(CMD);
% Load the fort.90 as well in misc folder
diff --git a/matlab/load_params.m b/matlab/load_params.m
index fe3e9925..7042b6c8 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -1,64 +1,63 @@
-CO = h5readatt(filename,'/data/input','CO');
-% K_N = h5readatt(filename,'/data/input','eta_n');
-% K_T = h5readatt(filename,'/data/input','eta_T');
-K_N = h5readatt(filename,'/data/input','K_n');
-K_T = h5readatt(filename,'/data/input','K_T');
-K_E = h5readatt(filename,'/data/input','K_E');
-Q0 = h5readatt(filename,'/data/input','q0');
-SHEAR = h5readatt(filename,'/data/input','shear');
-EPS = h5readatt(filename,'/data/input','eps');
-PMAXI = h5readatt(filename,'/data/input','pmaxi');
-JMAXI = h5readatt(filename,'/data/input','jmaxi');
-PMAXE = h5readatt(filename,'/data/input','pmaxe');
-JMAXE = h5readatt(filename,'/data/input','jmaxe');
-NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
-NU = h5readatt(filename,'/data/input','nu');
-Nx = h5readatt(filename,'/data/input','Nx');
-Ny = h5readatt(filename,'/data/input','Ny');
-L = h5readatt(filename,'/data/input','Lx');
-CLOS = h5readatt(filename,'/data/input','CLOS');
-DT_SIM = h5readatt(filename,'/data/input','dt');
-MU = h5readatt(filename,'/data/input','mu');
+DATA.CO = h5readatt(filename,'/data/input','CO');
+DATA.K_N = h5readatt(filename,'/data/input','eta_n');
+DATA.K_T = h5readatt(filename,'/data/input','eta_T');
+% DATA.K_N = h5readatt(filename,'/data/input','K_n');
+% DATA.K_T = h5readatt(filename,'/data/input','K_T');
+% DATA.K_E = h5readatt(filename,'/data/input','K_E');
+DATA.Q0 = h5readatt(filename,'/data/input','q0');
+DATA.SHEAR = h5readatt(filename,'/data/input','shear');
+DATA.EPS = h5readatt(filename,'/data/input','eps');
+DATA.PMAXI = h5readatt(filename,'/data/input','pmaxi');
+DATA.JMAXI = h5readatt(filename,'/data/input','jmaxi');
+DATA.PMAXE = h5readatt(filename,'/data/input','pmaxe');
+DATA.JMAXE = h5readatt(filename,'/data/input','jmaxe');
+DATA.NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
+DATA.NU = h5readatt(filename,'/data/input','nu');
+DATA.Nx = h5readatt(filename,'/data/input','Nx');
+DATA.Ny = h5readatt(filename,'/data/input','Ny');
+DATA.L = h5readatt(filename,'/data/input','Lx');
+DATA.CLOS = h5readatt(filename,'/data/input','CLOS');
+DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
+DATA.MU = h5readatt(filename,'/data/input','mu');
% MU = str2num(filename(end-18:end-14)); %bad...
-W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
-W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
-W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
-W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y';
-W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y';
+DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
+DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
+DATA.W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
+DATA.W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y';
+DATA.W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y';
-if NON_LIN == 'y'
- NON_LIN = 1;
+if DATA.NON_LIN == 'y'
+ DATA.NON_LIN = 1;
else
- NON_LIN = 0;
+ DATA.NON_LIN = 0;
end
-switch abs(CO)
- case 0; CONAME = 'LB';
- case 1; CONAME = 'DG';
- case 2; CONAME = 'SG';
- case 3; CONAME = 'PA';
- case 4; CONAME = 'FC';
- otherwise; CONAME ='UK';
+switch abs(DATA.CO)
+ case 0; DATA.CONAME = 'LB';
+ case 1; DATA.CONAME = 'DG';
+ case 2; DATA.CONAME = 'SG';
+ case 3; DATA.CONAME = 'PA';
+ case 4; DATA.CONAME = 'FC';
+ otherwise; DATA.CONAME ='UK';
end
-if (CO <= 0); CONAME = [CONAME,'DK'];
-else; CONAME = [CONAME,'GK'];
+if (DATA.CO <= 0); DATA.CONAME = [DATA.CONAME,'DK'];
+else; DATA.CONAME = [DATA.CONAME,'GK'];
end
-if (CLOS == 0); CLOSNAME = 'Trunc.';
-elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
-elseif(CLOS == 2); CLOSNAME = 'Clos. 2';
+if (DATA.CLOS == 0); DATA.CLOSNAME = 'Trunc.';
+elseif(DATA.CLOS == 1); DATA.CLOSNAME = 'Clos. 1';
+elseif(DATA.CLOS == 2); DATA.CLOSNAME = 'Clos. 2';
end
-if (PMAXE == PMAXI) && (JMAXE == JMAXI)
- degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)];
+if (DATA.PMAXE == DATA.PMAXI) && (DATA.JMAXE == DATA.JMAXI)
+ degngrad = ['P_',num2str(DATA.PMAXE),'_J_',num2str(DATA.JMAXE)];
else
- degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
- '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
+ degngrad = ['Pe_',num2str(DATA.PMAXE),'_Je_',num2str(DATA.JMAXE),...
+ '_Pi_',num2str(DATA.PMAXI),'_Ji_',num2str(DATA.JMAXI)];
end
degngrad = [degngrad,'_Kn_%1.1f_nu_%0.0e_',...
- CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[K_N,NU,MU]);
-if ~NON_LIN; degngrad = ['lin_',degngrad]; end
-resolution = [num2str(Nx),'x',num2str(Ny/2),'_'];
-gridname = ['L_',num2str(L),'_'];
-PARAMS = [resolution,gridname,degngrad];
-% BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
+ DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
+degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MU]);
+if ~DATA.NON_LIN; degngrad = ['lin_',degngrad]; end
+resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_'];
+gridname = ['L_',num2str(DATA.L),'_'];
+DATA.PARAMS = [resolution,gridname,degngrad];
\ No newline at end of file
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
new file mode 100644
index 00000000..f054b94d
--- /dev/null
+++ b/matlab/photomaton.m
@@ -0,0 +1,24 @@
+function [ FIGURE ] = photomaton( DATA,OPTIONS )
+%UNTITLED5 Summary of this function goes here
+%% Processing
+toplot = process_field(DATA,OPTIONS);
+FNAME = toplot.FILENAME;
+FRAMES = toplot.FRAMES;
+
+%
+TNAME = [];
+FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 numel(FRAMES) 1])
+ for i_ = 1:numel(FRAMES)
+ subplot(1,numel(FRAMES),i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))); set(pclr, 'edgecolor','none');pbaspect(toplot.ASPECT)
+ colormap(bluewhitered);
+ xlabel(toplot.XNAME); ylabel(toplot.YNAME);set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_))));
+ if OPTIONS.INTERP
+ shading interp;
+ end
+ end
+ legend(['$',toplot.FIELDNAME,'$']);
+ FIGURE.FIGNAME = [FNAME,'_snaps',TNAME];
+end
+
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index a8cdddeb..597952e0 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -1,73 +1,65 @@
+function [FIGURE] = plot_radial_transport_and_shear(DATA, TAVG_0, TAVG_1)
%Compute steady radial transport
tend = TAVG_1; tstart = TAVG_0;
-[~,its0D] = min(abs(Ts0D-tstart));
-[~,ite0D] = min(abs(Ts0D-tend));
-SCALE = (2*pi/Nx/Ny)^2;
-gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
-gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
-if W_HF
- q_infty_avg = mean(HFLUX_X(its0D:ite0D))*SCALE;
-end
-% Compute steady shearing rate
-tend = TAVG_1; tstart = TAVG_0;
-[~,its2D] = min(abs(Ts3D-tstart));
-[~,ite2D] = min(abs(Ts3D-tend));
-shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1));
-shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1));
-% plots
-fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
+[~,its0D] = min(abs(DATA.Ts0D-tstart));
+[~,ite0D] = min(abs(DATA.Ts0D-tend));
+SCALE = (2*pi/DATA.Nx/DATA.Ny)^2;
+gamma_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+gamma_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+
+FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
- yyaxis left
- plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
- plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
+% yyaxis left
+ plot(DATA.Ts0D,DATA.PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+ plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
- ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
- if W_HF
- yyaxis right
- plot(Ts0D,HFLUX_X*SCALE,'DisplayName','$\langle T_i \partial_y\phi \rangle_y$'); hold on;
- grid on; set(gca,'xticklabel',[]); ylabel('$Q_x$')
- xlim([Ts0D(1),Ts0D(end)]); ylim([0,5*abs(q_infty_avg)]);
+ ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+ title(['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ', $\kappa_N=$',num2str(DATA.K_N),...
+ ', $L=',num2str(DATA.L),'$, $N=',num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',num2str(DATA.JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(DATA.MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
+ %% zonal vs nonzonal energies for phi(t) and shear radial profile
+
+ % computation
+ Ns3D = numel(DATA.Ts3D);
+ [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ Ephi_Z = zeros(1,Ns3D);
+ Ephi_NZ_kgt0 = zeros(1,Ns3D);
+ Ephi_NZ_kgt1 = zeros(1,Ns3D);
+ Ephi_NZ_kgt2 = zeros(1,Ns3D);
+ high_k_phi = zeros(1,Ns3D);
+ dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+% dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ [amp,ikzf] = max(abs((KX~=0).*DATA.PHI(:,1,1,it)));
+ Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+ Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+ Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+ Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+ dxphi(:,:,1,it) = real(fftshift(ifft2(-1i*KX.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
+% dx2phi(:,:,1,it) = real(fftshift(ifft2(-KX.^2.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
end
- %
-% subplot(312)
-% clr = line_colors(1,:);
-% lstyle = line_styles(1);
-% % plt = @(x_) mean(x_,1);
-% plt = @(x_) x_(1,:);
-% plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(\partial^2_x\phi)$'); hold on;
-% plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle \partial^2_x\phi\rangle_y$'); hold on;
-% plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle \partial^2_x\phi\rangle_x$'); hold on;
-% plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle \partial^2_x\phi\rangle_{x,y}$'); hold on;
-% plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
-% 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
-% ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
-% grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
- %% zonal vs nonzonal energies for phi(t)
+ [~,its3D] = min(abs(DATA.Ts3D-tstart));
+ [~,ite3D] = min(abs(DATA.Ts3D-tend));
+ % plot
subplot(312)
it0 = 1; itend = Ns3D;
trange = it0:itend;
pltx = @(x) x;%-x(1);
- plty = @(x) x./max((x(1:end)));
+ plty = @(x) x./max((x(its3D:ite3D)));
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
- plot(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on;
- ylim([1e-1, 1.5]); ylabel('$E_{Z}$')
+ plot(pltx(DATA.Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',DATA.CONAME]); hold on;
+ ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
- plot(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$');
-% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName','$k_p>1$');
-% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName','$k_p>2$');
- % semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]);
-% legend('Location','southeast')
- xlim([Ts3D(it0), Ts3D(itend)]);
- ylim([1e-6, 1.0]); ylabel('$E_{NZ}$')
+ plot(pltx(DATA.Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$');
+ 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]);
subplot(313)
- [TY,TX] = meshgrid(x,Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
+ [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+ pclr = pcolor(TX,TY,squeeze((mean(dxphi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x\phi\rangle_y$') %colorbar;
+% pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
caxis([-3 3]);
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
-save_figure
+end
\ No newline at end of file
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index f30b11df..1d051ece 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -7,26 +7,6 @@ Ns3D = numel(Ts3D);
Ts5D = Ts5D';
Ts3D = Ts3D';
-%% Build grids
-Nkx = numel(kx); Nky = numel(ky);
-[KY,KX] = meshgrid(ky,kx);
-Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
-dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
-KPERP2 = KY.^2+KX.^2;
-[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
-[KY_XY,KX_XY] = meshgrid(ky,kx);
-[KZ_XZ,KX_XZ] = meshgrid(z,kx);
-[KZ_YZ,KY_YZ] = meshgrid(z,ky);
-
-Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z);
-Lx = 2*pi/dkx; Ly = 2*pi/dky;
-dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz;
-x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
-y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
-[Y_XY,X_XY] = meshgrid(y,x);
-[Z_XZ,X_XZ] = meshgrid(z,x);
-[Z_YZ,Y_YZ] = meshgrid(z,y);
-
%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Analysis :')
disp('- iFFT')
diff --git a/matlab/process_field.m b/matlab/process_field.m
new file mode 100644
index 00000000..563694e7
--- /dev/null
+++ b/matlab/process_field.m
@@ -0,0 +1,196 @@
+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
+%% 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);
+POLARPLOT = OPTIONS.POLARPLOT;
+LTXNAME = OPTIONS.NAME;
+switch OPTIONS.PLAN
+ case 'xy'
+ XNAME = '$x$'; YNAME = '$y$';
+ [Y,X] = meshgrid(DATA.y,DATA.x);
+ REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
+ case 'xz'
+ XNAME = '$x$'; YNAME = '$z$';
+ [Y,X] = meshgrid(DATA.z,DATA.x);
+ REALP = 1; COMPDIM = 2; SCALE = 0;
+ case 'yz'
+ XNAME = '$y$'; YNAME = '$z$';
+ [Y,X] = meshgrid(DATA.z,DATA.y);
+ REALP = 1; COMPDIM = 1; SCALE = 0;
+ case 'kxky'
+ XNAME = '$k_x$'; YNAME = '$k_y$';
+ [Y,X] = meshgrid(DATA.ky,DATA.kx);
+ 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 = 2; POLARPLOT = 0; SCALE = 0;
+ OPTIONS.COMP = 'max';
+ case 'kyz'
+ XNAME = '$k_y$'; YNAME = '$z$';
+ [Y,X] = meshgrid(DATA.z,DATA.ky);
+ REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
+ OPTIONS.COMP = 'max';
+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,:,:);
+ COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i));
+ FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
+ case 2
+ i = OPTIONS.COMP;
+ compr = @(x) x(:,i,:);
+ COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i));
+ FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
+ case 3
+ i = OPTIONS.COMP;
+ compr = @(x) x(:,:,i);
+ COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.z(i)/pi);
+ FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
+ end
+end
+
+switch REALP
+ case 1 % Real space plot
+ process = @(x) real(fftshift(ifft2(x,Nx,Ny)));
+ shift_x = @(x) x;
+ shift_y = @(x) x;
+ case 0 % Frequencies plot
+ OPTIONS.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 = FRAMES
+ tmp = squeeze(compr(DATA.PHI(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case 'n_i'
+ NAME = 'ni';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.DENS_I(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case 'n_i^{NZ}'
+ NAME = 'niNZ';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.DENS_I(:,:,:,it).*(KY~=0)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it).*(KY~=0)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+end
+TOPLOT.FIELD = FIELD;
+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;
+
+end
+
diff --git a/matlab/save_figure.m b/matlab/save_figure.m
index c5a9c72c..1aa1dd64 100644
--- a/matlab/save_figure.m
+++ b/matlab/save_figure.m
@@ -1,8 +1,10 @@
%% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME)
% and parameters
-if ~exist([BASIC.RESDIR,'/fig'], 'dir')
- mkdir([BASIC.RESDIR,'/fig'])
+function save_figure(data,FIGURE)
+if ~exist([data.localdir,'/fig'], 'dir')
+ mkdir([data.localdir,'/fig'])
end
-saveas(fig,[BASIC.RESDIR,'/fig/', FIGNAME,'.fig']);
-saveas(fig,[BASIC.RESDIR, FIGNAME,'.png']);
-disp(['Figure saved @ : ',[BASIC.RESDIR, FIGNAME,'.png']])
\ No newline at end of file
+saveas(FIGURE.fig,[data.localdir,'/fig/', FIGURE.FIGNAME,'.fig']);
+saveas(FIGURE.fig,[data.localdir, FIGURE.FIGNAME,'.png']);
+disp(['Figure saved @ : ',[data.localdir, FIGURE.FIGNAME,'.png']])
+end
\ No newline at end of file
diff --git a/matlab/show_geometry.m b/matlab/show_geometry.m
new file mode 100644
index 00000000..e57746d8
--- /dev/null
+++ b/matlab/show_geometry.m
@@ -0,0 +1,96 @@
+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);
+else % Zpinch geometry
+ Nturns = 0.1; Tgeom = 0;
+ a = 1; % 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, 3); % cylinder height
+ [t, p] = meshgrid(phi, theta);
+ x_tor = a.*cos(p);
+ z_tor = a.*sin(p);
+ y_tor = t;
+end
+FIGURE.fig = figure; set(gcf, 'Position', [50 50 1200 600])
+torus=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
+set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none')
+xlabel('X');ylabel('Y');zlabel('Z');
+% 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
+theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
+plot3(Xfl(theta),Yfl(theta),Zfl(theta))
+% Planes plot
+theta = DATA.z ; % Poloidal angle
+plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok')
+% Plot x,y,bvec at these points
+scale = 0.15;
+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');
+
+%% 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);
+max_ = 0;
+FIELDS = zeros(DATA.Nx,DATA.Ny,numel(OPTIONS.PLANES));
+
+for it_ = 1:numel(OPTIONS.TIME)
+[~,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);
+ 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
+end
+%%
+axis equal
+view([1,-2,1])
+
+end
+
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 5515dca9..249c2a42 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -2,33 +2,31 @@ addpath(genpath('../matlab')) % ... add
addpath(genpath('../matlab/plots')) % ... add
outfile ='';
%% Directory of the simulation
-if 1% Local results
+if 0% Local results
outfile ='';
-outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_adiabe';
-% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_lin';
-% outfile ='fluxtube_salphaB_s0/64x64x16_5x3_L_200_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK';
-% outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
-% outfile ='Linear_Device/64x64x20_5x2_Lx_20_Ly_150_q0_25_kN_0.24_kT_0.03_nu_1e-02_DGGK';
- BASIC.RESDIR = ['../results/',outfile,'/'];
- BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
- system(['mkdir -p ',BASIC.MISCDIR]);
- CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='fluxtube_salphaB_s0/100x100x30_5x3_Lx_200_Ly_100_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe';
+% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe';
+% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe_Sg';
+ LOCALDIR = ['../results/',outfile,'/'];
+ MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
+ system(['mkdir -p ',MISCDIR]);
+ CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x200_L_200_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
-BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
-BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
+% BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
+MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
-% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A
-compile_results %Compile the results from first output found to JOBNUMMAX if existing
+JOBNUMMIN = 00; JOBNUMMAX = 20;
+data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-%% Post-processing
-post_processing
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
@@ -37,155 +35,48 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1500; TAVG_1 = 4000; % Averaging times duration
-plot_radial_transport_and_shear
+TAVG_0 = 1500; TAVG_1 = 5000; % Averaging times duration
+fig = plot_radial_transport_and_shear(data,TAVG_0,TAVG_1);
+save_figure(data,fig)
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-t0 =000; iz = 1; ix = 1; iy = 1;
-skip_ =1; FPS = 30; DELAY = 1/FPS;
-[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
-[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
-T = Ts3D; FRAMES = FRAMES_3D;
-INTERP = 0; NORMALIZED = 1; CONST_CMAP = 0; BWR =1;% Gif options
-% Field to plot
-% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i';
-% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
-% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
-% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
-% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
-% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
-% FIELD = phi-Z_phi; NAME = 'NZphi'; FIELDNAME = '\phi_{NZ}';
-% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
-
-% Sliced
-% plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
-% plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
-plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
-
-% K-space
-% FIELD = PHI.*(KY~=0); NAME = 'NZPHI'; FIELDNAME = '\tilde \phi_{k_y\neq0}';
-% FIELD = Ne00; NAME = 'Ne00'; FIELDNAME = 'N_e^{00}';
-% FIELD = Ni00; NAME = 'Ni00'; FIELDNAME = 'N_i^{00}';
-% plt = @(x) fftshift((abs(x( :, :,1,:))),2); X = fftshift(KX,2); Y = fftshift(KY,2); XNAME = 'k_x'; YNAME = 'k_y';
-
-% Averaged
-% plt = @(x) mean(x,1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
-% plt = @(x) mean(x,2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
-% plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
-
-FIELD = squeeze(plt(FIELD));
-
-% Naming
-GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS];
-
-% Create movie (gif or mp4)
-create_gif
-% create_mov
+options.INTERP = 1;
+options.POLARPLOT = 0;
+options.NAME = 'n_i^{NZ}';
+options.PLAN = 'xy';
+options.COMP = 1;
+options.TIME = 00:20:6000;
+% options.TIME = 140:0.5:160;
+data.a = data.EPS * 2000;
+create_film(data,options,'.gif')
end
if 0
-%% 2D plot : real space
-
-% Chose the field to plot
-% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
-% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
-% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
-% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
-% FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}';
-% FIELD = dens_i-Z_n_i; FNAME = 'ni_NZ'; FIELDLTX = 'n_i^{NZ}';
-% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
-% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
-% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
-FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
-% FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}';
-% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
-% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
-
-% Chose when to plot it
-tf = [0 15 27 28 30];
-
-% Planar plot: choose a plane to plot at x0/y0/z0 coordinates
-x0 = 0.0; y0 = 0.0; z0 = 0.0;
-[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0));
-% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')']
-% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')']
-plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(z=',num2str(round(z(iz)/pi)),'\pi)']
-
-% Averaged
-% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_x']
-% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_y']
-% plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_z']
-
-
-%
-TNAME = [];
-fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
-plt_2 = @(x) x;%./max(max(x));
- for i_ = 1:numel(tf)
- [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
- subplot(1,numel(tf),i_)
- DATA = plt_2(squeeze(plt(FIELD,it)));
- pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap(bluewhitered); %caxis([-30,30]);
- xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
- end
- legend(latexize(FIELDLTX));
-save_figure
+%% 2D snapshots
+% Options
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.NAME = 'n_i^{NZ}';
+options.PLAN = 'xy';
+options.COMP = 1;
+options.TIME = [100 400 2500 5500];
+data.a = data.EPS * 1000;
+fig = photomaton(data,options);
+save_figure(data,fig)
end
if 0
-%% 2D plot : k space
-
-% Chose the field to plot
-% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
-% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
-
-% Chose when to plot it
-% tf = [0 15 27 28 30];
-tf = 200:200:1200;
-% tf = 8000;
-
-% Planar plot: choose a plane to plot at x0/y0/z0 coordinates
-x0 = 0.0; y0 = 0.3; z0 = 0.5*pi;
-[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0));
-% plt = @(x,it) abs(x(ix, :, :,it)); X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_x=',num2str(round(kx(ix))),')'];
-% plt = @(x,it) abs(x( :,iy, :,it)); X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_y=',num2str(round(ky(iy))),')'];
-% plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)'];
-% %
-% plt_x = @(x) fftshift(x,1); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,1); plt = @(x,it) max(abs(x( :, :,:,it)),[],1);
-% X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_x)'];
-%
-% plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],2);
-% X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_y)'];
-
-plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,2); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],3);
-X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(\max_z)'];
-
-%
-TNAME = [];
-fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 300*numel(tf), 500])
- for i_ = 1:numel(tf)
- [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
- subplot(1,numel(tf),i_)
- DATA = plt_2(squeeze(plt(FIELD,it)));
- pclr = pcolor(plt_x(X),plt_y(Y),plt_z(DATA)); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
- caxis([0 1]*1e3);
-% caxis([-1 1]*5);
- xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
- if(i_ > 1); set(gca,'ytick',[]); end;
- title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
- end
- legend(latexize(FIELDLTX));
-save_figure
+%% 3D plot on the geometry
+options.INTERP = 1;
+options.NAME = 'n_i';
+options.PLANES = 1;
+options.TIME = 5000;
+data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio
+FIGURE = show_geometry(data,options);
end
if 0
@@ -319,35 +210,4 @@ FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0;
plot_kperp_spectrum
-end
-
-if 0
-%% Torus plot
-aminor = EPS; % Torus minor radius
-Rmajor = 1.; % Torus major radius
-theta = linspace(-pi, pi, 30) ; % Poloidal angle
-phi = linspace(0., 2.*pi, 30) ; % Toroidal angle
-[t, p] = meshgrid(phi, theta);
-x = (Rmajor + aminor.*cos(p)) .* cos(t);
-y = (Rmajor + aminor.*cos(p)) .* sin(t);
-z = aminor.*sin(p);
-figure;
-torus=surf(x, y, z); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
-set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none')
-xlabel('X');ylabel('Y');zlabel('Z');
-% field line plot
-Nturns = 1;
-theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
-xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0);
-yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0);
-zFL = aminor.*sin(theta);
-plot3(xFL,yFL,zFL,'r')
-% Planes plot
-theta = linspace(-pi, pi, Nz) ; % Poloidal angle
-xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0);
-yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0);
-zFL = aminor.*sin(theta);
-plot3(xFL,yFL,zFL,'sb')
-axis equal
-
end
\ No newline at end of file
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index 95365b4d..7f0d2c3e 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -13,16 +13,17 @@ K_T = 0.0; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-NX = 100; % real space x-gridpoints
+NX = 150; % real space x-gridpoints
NY = 1; % '' y-gridpoints
-LX = 150; % Size of the squared frequency domain
+LX = 200; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
+SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
+TMAX = 200; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
@@ -36,7 +37,7 @@ NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
-CO = 2;
+CO = 4;
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
@@ -89,9 +90,9 @@ for i = 1:Nparam
system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
% Load and process results
%%
diff --git a/wk/local_run.m b/wk/local_run.m
index 10fefe6e..5490687b 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -12,11 +12,11 @@ SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NU_HYP = 0.0;
KIN_E = 0; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
-NX = 50; % Spatial radial resolution ( = 2x radial modes)
-LX = 300; % Radial window size
+NX = 100; % Spatial radial resolution ( = 2x radial modes)
+LX = 200; % Radial window size
NY = 100; % Spatial azimuthal resolution (= azim modes)
-LY = 300; % Azimuthal window size
-NZ = 20; % number of perpendicular planes (parallel grid)
+LY = 100; % Azimuthal window size
+NZ = 30; % number of perpendicular planes (parallel grid)
P = 4;
J = 2;
%% GEOMETRY PARAMETERS
@@ -27,8 +27,8 @@ GRADB = 1.0; % Magnetic gradient
CURVB = 1.0; % Magnetic curvature
SG = 1; % Staggered z grids option
%% TIME PARAMETERS
-TMAX = 50; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 200; % Maximal time unit
+DT = 5e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 3D arrays
@@ -40,7 +40,7 @@ JOB2LOAD= -1;
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
CO = 1;
CLOS = 0; % Closure model (0: =0 truncation)
-NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
SIMID = 'fluxtube_salphaB_s0'; % Name of the simulation
% SIMID = 'simulation_A'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 15674da9..e6ea6a5e 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -3,8 +3,7 @@ addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain)
% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_3.9';
-for K_N = [1/0.6]
+ EXECNAME = 'helaz_3.5';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -20,33 +19,34 @@ NP_KX = 24; % MPI processes along kx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
-K_N = 1.0/0.6; % Density gradient drive (R/Ln)
+K_N = 1.0/0.6; % Density gradient drive (R/Ln)
NU_HYP = 0.0;
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-Nx = 300; % Realspace x-gridpoints
-Ny = 300; % Realspace y-gridpoints
-Lx = 120; % Size of the squared frequency domain
-Ly = 120; % Size of the squared frequency domain
-Nz = 1; % number of perpendicular planes (parallel grid)
-q0 = 1.0; % q factor ()
-shear = 0.0; % magnetic shear
-eps = 0.0; % inverse aspect ratio
-P = 8;
-J = 4;
+NX = 300; % Realspace x-gridpoints
+NY = 300; % Realspace y-gridpoints
+LX = 120; % Size of the squared frequency domain
+LY = 120; % Size of the squared frequency domain
+NZ = 1; % number of perpendicular planes (parallel grid)
+Q0 = 1.0; % q factor ()
+SHEAR = 0.0; % magnetic shear
+EPS = 0.0; % inverse aspect ratio
+P = 4;
+J = 2;
%% TIME PARAMETERS
TMAX = 10000; % Maximal time unit
-DT = 8e-3; % Time step
+DT = 2e-3; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS3D = 1/2; % Sampling per time unit for 3D arrays
-SPS5D = 1/100; % Sampling per time unit for 5D arrays
-JOB2LOAD= 2; % start from t=0 if <0, else restart from outputs_$job2load
+SPS5D = 1/50; % Sampling per time unit for 5D arrays
+JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO = 3;
+CO = 2;
CLOS = 0; % Closure model (0: =0 truncation)
-NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
% SIMID = 'test_chained_job'; % Name of the simulation
SIMID = 'simulation_A'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
@@ -63,6 +63,10 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% unused
+KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
+GRADB = 1.0; % Magnetic gradient
+CURVB = 1.0; % Magnetic curvature
+SG = 0; % Staggered z grids option
PMAXE = P; % Highest electron Hermite polynomial degree
JMAXE = J; % Highest '' Laguerre ''
PMAXI = P; % Highest ion Hermite polynomial degree
@@ -72,14 +76,15 @@ KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KXEQ0 = 0; % put kx = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = Nx*pi/Lx;% Highest fourier mode
+kmax = NX*pi/LX;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
% kmaxcut = 2.5;
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
BCKGD0 = 0.0; % Init background
TAU = 1.0; % e/i temperature ratio
-K_T = 0.0; % Temperature gradient
+K_T = 0.0; % Temperature gradient
+K_E = 0.0; % ES '''
INIT_PHI= 1; % Start simulation with a noisy phi and moments
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
@@ -115,4 +120,3 @@ if(SUBMIT)
end
system('rm fort*.90');
disp('done');
-end
diff --git a/wk/photomaton.m b/wk/photomaton.m
deleted file mode 100644
index 39f6813d..00000000
--- a/wk/photomaton.m
+++ /dev/null
@@ -1,59 +0,0 @@
-if 0
-%% Photomaton : real space
-
-% Chose the field to plot
-% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
-% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
-% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
-% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
-% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
-% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
-FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
-
-% Chose when to plot it
-tf = [0 10 50 100];
-
-% Slice
-ix = 1; iy = 1; iz = 1;
-% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')']
-% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')']
-% plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(z(iz))),')']
-
-% Averaged
-% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_x']
-% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_y']
-plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_z']
-
-
-%
-TNAME = [];
-fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
-plt_2 = @(x) x./max(max(x));
- for i_ = 1:numel(tf)
- [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
- subplot(1,numel(tf),i_)
- DATA = plt_2(squeeze(plt(FIELD,it)));
- pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap(bluewhitered); caxis([-1,1]);
- xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
- end
- legend(latexize(FIELDLTX));
-save_figure
-end
-
-if 0
-%% Photomaton : quiver ExB velocity
-figure
-skip = 2;
-plt = @(x) x./max(max(x));
-FNAME = 'ZF'; FIELDLTX = '$\bm{u}^{ZF}$'; X_XY = RR; Y_XY = ZZ;
-tf = 1200; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-UY = plt(drphi(1:skip:end,1:skip:end,it1)); UX = plt(-dzphi(1:skip:end,1:skip:end,it1));
-pclr = pcolor(X_XY,Y_XY,plt(ni00(:,:,it1))); set(pclr, 'edgecolor','none');
-hold on
-quiver((X_XY(1:skip:end,1:skip:end)),(Y_XY(1:skip:end,1:skip:end)),UX,UY,'r'); xlim(L/2*[-1 1]); ylim(L/2*[-1 1]);
-pbaspect([1 1 1])
-xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
-title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
-end
\ No newline at end of file
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
index 87004414..a86a5247 100644
--- a/wk/shearless_linear_fluxtube.m
+++ b/wk/shearless_linear_fluxtube.m
@@ -13,17 +13,18 @@ K_T = 6.0; % Temperature '''
SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % Kinetic (1) or adiabatic (0) electron model
%% GRID PARAMETERS
-NX = 1; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
-LX = 0; % Size of the squared frequency domain
+NX = 2; % real space x-gridpoints
+NY = 16; % '' y-gridpoints
+LX = 2*pi/0.25; % Size of the squared frequency domain
LY = 2*pi/0.25; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
Q0 = 2.7; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.18; % inverse aspect ratio
+SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 10; % Maximal time unit
-DT = 1e-3; % Time step
+TMAX = 40; % Maximal time unit
+DT = 8e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 10; % Sampling per time unit for 2D arrays
@@ -49,7 +50,7 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-kmax = Nx*pi/Lx;% Highest fourier mode
+kmax = NX*pi/LX;% Highest fourier mode
NU_HYP = 0.0; % Hyperdiffusivity coefficient
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
@@ -65,7 +66,6 @@ KXEQ0 = 0; % put kx = 0
NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
%% PARAMETER SCANS
-if 1
%% Parameter scan over PJ
% PA = [2 4];
% JA = [1 2];
@@ -77,11 +77,11 @@ mup_ = MU_P;
muj_ = MU_J;
Nparam = numel(PA);
param_name = 'PJ';
-gamma_Ni00 = zeros(Nparam,floor(Nx/2)+1);
-gamma_Nipj = zeros(Nparam,floor(Nx/2)+1);
-gamma_phi = zeros(Nparam,floor(Nx/2)+1);
-% Ni00_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D);
-% PHI_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D);
+Nkx = NX;
+Nky = NY/2+1;
+gamma_Ni00 = zeros(Nparam,Nkx,Nky);
+gamma_Nipj = zeros(Nparam,Nkx,Nky);
+gamma_phi = zeros(Nparam,Nkx,Nky);
for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
@@ -90,21 +90,68 @@ for i = 1:Nparam
setup
% system(['rm fort*.90']);
% Run linear simulation
- if RUN
+ if 0
system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
% Load and process results
%%
filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
load_results
+
+ for ikx = 1:Nkx
+ for iky = 1:Nky
+ tend = max(Ts3D(abs(Ni00(ikx,iky,1,:))~=0));
+ tstart = 0.2*tend;
+ [~,itstart] = min(abs(Ts3D-tstart));
+ [~,itend] = min(abs(Ts3D-tend));
+ trange = itstart:itend;
+ % exp fit on moment 00
+ X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,iky,1,trange)));
+ gamma_Ni00(i,ikx,iky) = LinearFit_s(X_,Y_);
+ % exp fit on phi
+ X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,iky,1,trange)));
+ gamma_phi (i,ikx,iky) = LinearFit_s(X_,Y_);
+ end
+ end
+ gamma_Ni00(i,:,:) = real(gamma_Ni00(i,:,:));% .* (gamma_Ni00(i,:)>=0.0));
+ gamma_Nipj(i,:,:) = real(gamma_Nipj(i,:,:));% .* (gamma_Nipj(i,:)>=0.0));
+ if 0
+ %% Fit verification
+ figure;
+ for i = 1:1:Nky
+ X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(1,i,1,:)));
+ semilogy(X_,Y_,'DisplayName',['k=',num2str(ky(i))]); hold on;
+ end
+ end
end
+if 1
+%% Plot
+SCALE = 1;%sqrt(2);
+fig = figure; FIGNAME = 'linear_study';
+plt = @(x) squeeze(x);
+% subplot(211)
+ for i = 1:Nparam
+ clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
+ linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
+ plot(plt(SCALE*ky(1:Nky)),plt(gamma_phi(i,2,1:end)),...
+ 'Color',clr,...
+ 'LineStyle',linestyle{1},'Marker','^',...
+...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
+ 'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']);
+ hold on;
+ end
+ grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(ky)]);
+ title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
+ legend('show'); %xlim([0.01,10])
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']);
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']);
end
if 0
%% Trajectories of some modes
figure;
-for i = 1:10:Nx/2+1
- semilogy(Ts3D,squeeze(abs(Ne00(i,2,1,:))),'DisplayName',['k=',num2str(kx(i))]); hold on;
+for i = 1:1:Nky
+ semilogy(Ts3D,squeeze(abs(Ni00(1,i,1,:))),'DisplayName',['k=',num2str(ky(i))]); hold on;
end
end
\ No newline at end of file
--
GitLab