Skip to content
Snippets Groups Projects
Commit ca8720fd authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

faster and more versatile matlab scripts

parent 6d84fe59
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 286 additions and 379 deletions
matlab/load/compile_results_3Da.m 0 → 100644
+ 37
− 0
View file @ ca8720fd
function [field, Ts3D] = compile_results_3Da(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname)
CONTINUE = 1;
JOBNUM = JOBNUMMIN; JOBFOUND = 0;
% field
field = [];
Ts3D = [];
ii = 1;
while(CONTINUE)
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
% Check presence and jobnummax
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
%test if it is corrupted or currently running
try
openable = ~isempty(h5read(filename,'/data/var3d/time'));
catch
openable = 0;
end
if openable
% load field %%
[ F, T, ~] = load_3Da_data(filename, fieldname);
field = cat(5,field,F);
Ts3D = cat(1,Ts3D,T);
ii = ii + 1;
JOBFOUND = JOBFOUND + 1;
end
elseif (JOBNUM > JOBNUMMAX)
CONTINUE = 0;
end
JOBNUM = JOBNUM + 1;
end
if(JOBFOUND == 0)
disp('no results found, please verify the paths');
return;
end
end
\ No newline at end of file
matlab/load/compile_results_low_mem.m
+ 48
− 47
View file @ ca8720fd
function [DATA] = compile_results_low_mem(DATA,DIRECTORY,JOBNUMMIN,JOBNUMMAX)
CONTINUE = 1;
JOBNUM = JOBNUMMIN; JOBFOUND = 0;
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.MUx_EVOL = []; % evolution of parameter mu between jobs
DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
DATA.MUz_EVOL = []; % evolution of parameter mu between jobs
DATA.K_N_EVOL = []; %
DATA.K_T_EVOL = []; %
DATA.L_EVOL = []; %
DATA.DT_EVOL = []; %
DATA.params_evol.TJOB_SE = []; % Start and end times of jobs
DATA.params_evol.NU_EVOL = []; % evolution of parameter nu between jobs
DATA.params_evol.CO_EVOL = []; % evolution of CO
DATA.params_evol.MUx_EVOL = []; % evolution of parameter mu between jobs
DATA.params_evol.MUy_EVOL = []; % evolution of parameter mu between jobs
DATA.params_evol.MUz_EVOL = []; % evolution of parameter mu between jobs
DATA.params_evol.K_N_EVOL = []; %
DATA.params_evol.K_T_EVOL = []; %
DATA.params_evol.L_EVOL = []; %
DATA.params_evol.DT_EVOL = []; %
% Low memoery cost data
HFLUXI_ = [];
HFLUX_ = [];
HFLUXE_ = [];
GGAMMAI_ = [];
PGAMMAI_ = [];
GGAMMA_ = [];
PGAMMA_ = [];
Ts0D_ = [];
DATA.outfilenames = {};
ii = 1;
......@@ -36,38 +36,38 @@ while(CONTINUE)
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
[Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
[Parray, Jarray, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input/diag_par','write_gamma'),'y');
W_HF = strcmp(h5readatt(filename,'/data/input/diag_par','write_hf' ),'y');
KIN_E = strcmp(h5readatt(filename,'/data/input/model', 'ADIAB_E' ),'n');
BETA = h5readatt(filename,'/data/input/model','beta');
if W_GAMMA
[ GGAMMA_RI, Ts0D, ~] = load_0D_data(filename, 'gflux_xi');
PGAMMA_RI = load_0D_data(filename, 'pflux_xi');
GGAMMAI_ = cat(1,GGAMMAI_,GGAMMA_RI); clear GGAMMA_RI
PGAMMAI_ = cat(1,PGAMMAI_,PGAMMA_RI); clear PGAMMA_RI
[ GGAMMA_R, Ts0D, ~] = load_0D_data(filename, 'gflux_x');
PGAMMA_R = load_0D_data(filename, 'pflux_x');
GGAMMA_ = cat(1,GGAMMA_,GGAMMA_R); clear GGAMMA_R
PGAMMA_ = cat(1,PGAMMA_,PGAMMA_R); clear PGAMMA_R
end
if W_HF
[ HFLUX_XI, Ts0D, ~] = load_0D_data(filename, 'hflux_xi');
HFLUXI_ = cat(1,HFLUXI_,HFLUX_XI); clear HFLUX_XI
[ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x');
HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X
end
Ts0D_ = cat(1,Ts0D_,Ts0D);
% Evolution of simulation parameters
load_params
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.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
DATA.MUz_EVOL = [DATA.MUz_EVOL DATA.MUz DATA.MUz];
DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
DATA.K_T_EVOL = [DATA.K_T_EVOL DATA.K_T DATA.K_T];
DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
DATA = load_params(DATA,filename);
DATA.params_evol.TJOB_SE = [DATA.params_evol.TJOB_SE Ts0D(1) Ts0D(end)];
DATA.params_evol.NU_EVOL = [DATA.params_evol.NU_EVOL DATA.inputs.NU DATA.inputs.NU];
DATA.params_evol.CO_EVOL = [DATA.params_evol.CO_EVOL DATA.inputs.CO DATA.inputs.CO];
DATA.params_evol.MUx_EVOL = [DATA.params_evol.MUx_EVOL DATA.inputs.MUx DATA.inputs.MUx];
DATA.params_evol.MUy_EVOL = [DATA.params_evol.MUy_EVOL DATA.inputs.MUy DATA.inputs.MUy];
DATA.params_evol.MUz_EVOL = [DATA.params_evol.MUz_EVOL DATA.inputs.MUz DATA.inputs.MUz];
DATA.params_evol.K_N_EVOL = [DATA.params_evol.K_N_EVOL DATA.inputs.K_N DATA.inputs.K_N];
DATA.params_evol.K_T_EVOL = [DATA.params_evol.K_T_EVOL DATA.inputs.K_T DATA.inputs.K_T];
DATA.params_evol.L_EVOL = [DATA.params_evol.L_EVOL DATA.inputs.L DATA.inputs.L];
DATA.params_evol.DT_EVOL = [DATA.params_evol.DT_EVOL DATA.inputs.DT_SIM DATA.inputs.DT_SIM];
ii = ii + 1;
JOBFOUND = JOBFOUND + 1;
......@@ -118,30 +118,31 @@ else
% scaling
DATA.scale = 1;%(1/Nx/Ny)^2;
% Fields
DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_;
if(KIN_E)
DATA.HFLUX_XE = HFLUXE_;
end
DATA.Ts0D = Ts0D_;
DATA.KIN_E=KIN_E;
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
DATA.Je = Je; DATA.Ji = Ji;
DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.Npol = -z(1)/pi;
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.grids.Parray = Parray;
DATA.grids.Jarray = Jarray;
DATA.grids.kx = kx; DATA.grids.ky = ky; DATA.grids.z = z; DATA.grids.Npol = -z(1)/pi;
DATA.grids.x = x; DATA.grids.y = y;
DATA.grids.ikx0 = ikx0; DATA.grids.iky0 = iky0;
DATA.grids.Nx = Nx; DATA.grids.Ny = Ny; DATA.grids.Nz = Nz; DATA.grids.Nkx = Nkx; DATA.grids.Nky = Nky;
DATA.grids.Np = numel(Parray);DATA.grids.Nj = numel(Jarray);
DATA.grids.deltap = Parray(2)-Parray(1);
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
', $\kappa_{Ni}=$',num2str(DATA.K_N),', $\kappa_{Ti}=$',num2str(DATA.K_T),...
', $L=',num2str(DATA.L),'$, $N=',...
num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.Pi(end)),',',...
num2str(DATA.Ji(end)),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
',',num2str(DATA.MUy),')'];
DATA.paramshort = [num2str(DATA.Pmaxi),'x',num2str(DATA.Jmaxi),'x',...
num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz)];
DATA.param_title=['$\nu_{',DATA.inputs.CONAME,'}=$', num2str(DATA.inputs.NU), ...
', $\kappa_{Ni}=$',num2str(DATA.inputs.K_N),', $\kappa_{Ti}=$',num2str(DATA.inputs.K_T),...
', $L=',num2str(DATA.inputs.L),'$, $N=',...
num2str(DATA.grids.Nx),'$, $(P,J)=(',num2str(DATA.grids.Parray(end)),',',...
num2str(DATA.grids.Jarray(end)),')$,',' $\mu_{hd}=$(',num2str(DATA.inputs.MUx),...
',',num2str(DATA.inputs.MUy),')'];
DATA.paramshort = [num2str(DATA.grids.Np),'x',num2str(DATA.grids.Nj),'x',...
num2str(DATA.grids.Nkx),'x',num2str(DATA.grids.Nky),'x',num2str(DATA.grids.Nz)];
JOBNUM = LASTJOB;
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
......
matlab/load/load_3D_data.m
+ 16
− 6
View file @ ca8720fd
......@@ -13,6 +13,10 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
sz = size(tmp);
cmpx = 0;
end
% add a z dimension even if 2D
if(numel(sz) == 2)
sz = [sz 1];
end
% add time dimension
sz = [sz numel(time)];
data = zeros(sz);
......@@ -20,16 +24,22 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
if cmpx
if(numel(sz) == 3)
data(:,:,it) = tmp.real + 1i * tmp.imaginary;
else
switch numel(sz)
case(3)
data(:,:,1,it) = tmp.real + 1i * tmp.imaginary;
case(4)
data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
case(5)
data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end
else
if(numel(sz) == 3)
data(:,:,it) = tmp;
else
switch numel(sz)
case(3)
data(:,:,1,it) = tmp;
case(4)
data(:,:,:,it) = tmp;
case(5)
data(:,:,:,:,it) = tmp;
end
end
end
......
matlab/load/load_3Da_data.m 0 → 100644
+ 43
− 0
View file @ ca8720fd
function [ data, time, dt ] = load_3Da_data( filename, variablename )
%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
time = h5read(filename,'/data/var3d/time');
dt = h5readatt(filename,'/data/input/basic','dt');
cstart= h5readatt(filename,'/data/input/basic','start_iframe3d');
Na = h5readatt(filename,'/data/input/model','Na');
Np = h5readatt(filename,'/data/input/grid', 'Np');
Nj = h5readatt(filename,'/data/input/grid', 'Nj');
Nky = h5readatt(filename,'/data/input/grid', 'Nky');
Nkx = h5readatt(filename,'/data/input/grid', 'Nkx');
Nz = h5readatt(filename,'/data/input/grid', 'Nz');
% Find array size by loading the first output
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+1,'%06d')]);
try % check if it is complex or real
sz = size(tmp.real);
cmpx = 1;
catch
sz = size(tmp);
cmpx = 0;
end
% add a dimension if nz=1 or na=1
% if Na == 1
% sz = [1 sz];
% end
if Nz == 1
sz = [sz 1];
end
% add time dimension
data = zeros([sz numel(time)]);
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
if cmpx
data(:,:,:,:,it) = reshape(tmp.real,sz) + 1i * reshape(tmp.imaginary,sz);
else
data(:,:,:,:,it) = reshape(tmp,sz);
end
end
end
matlab/load/load_grid_data.m
+ 3
− 5
View file @ ca8720fd
function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data( filename )
function [ p, j, kx, ky, z ] = load_grid_data( filename )
%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ
pe = h5read(filename,'/data/grid/coordp');
je = h5read(filename,'/data/grid/coordj');
pi = h5read(filename,'/data/grid/coordp');
ji = h5read(filename,'/data/grid/coordj');
p = h5read(filename,'/data/grid/coordp');
j = h5read(filename,'/data/grid/coordj');
kx = h5read(filename,'/data/grid/coordkx');
ky = h5read(filename,'/data/grid/coordky');
z = h5read(filename,'/data/grid/coordz');
......
matlab/load/load_params.m
+ 75
− 73
View file @ ca8720fd
DATA.CO = h5readatt(filename,'/data/input/coll','CO');
DATA.K_N = h5readatt(filename,'/data/input/ions','k_N');
DATA.K_T = h5readatt(filename,'/data/input/ions','k_T');
DATA.Q0 = h5readatt(filename,'/data/input/geometry','q0');
DATA.EPS = h5readatt(filename,'/data/input/geometry','eps');
DATA.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
DATA.GEOM = h5readatt(filename,'/data/input/geometry','geometry');
% DATA.KAPPA = h5readatt(filename,'/data/input/geometry','kappa');
% DATA.DELTA = h5readatt(filename,'/data/input/geometry','delta');
function [DATA] = load_params(DATA,filename)
DATA.inputs.CO = h5readatt(filename,'/data/input/coll','CO');
DATA.inputs.K_N = h5readatt(filename,'/data/input/ions','k_N');
DATA.inputs.K_T = h5readatt(filename,'/data/input/ions','k_T');
DATA.inputs.Q0 = h5readatt(filename,'/data/input/geometry','q0');
DATA.inputs.EPS = h5readatt(filename,'/data/input/geometry','eps');
DATA.inputs.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
DATA.inputs.GEOM = h5readatt(filename,'/data/input/geometry','geometry');
% DATA.KAPPA = h5readatt(filename,'/data/input/geometry','kappa');
% DATA.DELTA = h5readatt(filename,'/data/input/geometry','delta');
DATA.DT_SIM = h5readatt(filename,'/data/input/basic','dt');
DATA.PMAX = h5readatt(filename,'/data/input/grid','pmax');
DATA.JMAX = h5readatt(filename,'/data/input/grid','jmax');
DATA.Nx = h5readatt(filename,'/data/input/grid','Nx');
DATA.Ny = h5readatt(filename,'/data/input/grid','Ny');
DATA.L = h5readatt(filename,'/data/input/grid','Lx');
try
DATA.CLOS = h5readatt(filename,'/data/input/model','CLOS');
DATA.NL_CLOS = h5readatt(filename,'/data/input/model','NL_CLOS');
catch
DATA.inputs.DT_SIM = h5readatt(filename,'/data/input/basic','dt');
DATA.inputs.PMAX = h5readatt(filename,'/data/input/grid','pmax');
DATA.inputs.JMAX = h5readatt(filename,'/data/input/grid','jmax');
DATA.inputs.Nx = h5readatt(filename,'/data/input/grid','Nx');
DATA.inputs.Ny = h5readatt(filename,'/data/input/grid','Ny');
DATA.inputs.L = h5readatt(filename,'/data/input/grid','Lx');
try
DATA.ha_cl = h5readatt(filename,'/data/input/closure','hierarchy_closure');
DATA.CLOS = h5readatt(filename,'/data/input/closure','dmax');
DATA.nl_cl = h5readatt(filename,'/data/input/closure','nonlinear_closure');
DATA.NL_CLOS = h5readatt(filename,'/data/input/closure','nmax');
DATA.inputs.CLOS = h5readatt(filename,'/data/input/model','CLOS');
DATA.inputs.NL_CLOS = h5readatt(filename,'/data/input/model','NL_CLOS');
catch
DATA.CLOS = 99;
DATA.NL_CLOS = 99;
try
DATA.inputs.ha_cl = h5readatt(filename,'/data/input/closure','hierarchy_closure');
DATA.inputs.CLOS = h5readatt(filename,'/data/input/closure','dmax');
DATA.inputs.nl_cl = h5readatt(filename,'/data/input/closure','nonlinear_closure');
DATA.inputs.NL_CLOS = h5readatt(filename,'/data/input/closure','nmax');
catch
DATA.inputs.CLOS = 99;
DATA.inputs.NL_CLOS = 99;
end
end
end
DATA.Na = h5readatt(filename,'/data/input/model','Na');
DATA.NU = h5readatt(filename,'/data/input/model','nu');
DATA.MUp = h5readatt(filename,'/data/input/model','mu_p');
DATA.MUj = h5readatt(filename,'/data/input/model','mu_j');
DATA.MUx = h5readatt(filename,'/data/input/model','mu_x');
DATA.MUy = h5readatt(filename,'/data/input/model','mu_y');
DATA.MUz = h5readatt(filename,'/data/input/model','mu_z');
DATA.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
DATA.BETA = h5readatt(filename,'/data/input/model','beta');
DATA.TAU_E = h5readatt(filename,'/data/input/model','tau_e');
DATA.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
DATA.K_cB = h5readatt(filename,'/data/input/model','k_cB');
DATA.K_gB = h5readatt(filename,'/data/input/model','k_gB');
DATA.inputs.Na = h5readatt(filename,'/data/input/model','Na');
DATA.inputs.NU = h5readatt(filename,'/data/input/model','nu');
DATA.inputs.MUp = h5readatt(filename,'/data/input/model','mu_p');
DATA.inputs.MUj = h5readatt(filename,'/data/input/model','mu_j');
DATA.inputs.MUx = h5readatt(filename,'/data/input/model','mu_x');
DATA.inputs.MUy = h5readatt(filename,'/data/input/model','mu_y');
DATA.inputs.MUz = h5readatt(filename,'/data/input/model','mu_z');
DATA.inputs.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
DATA.inputs.BETA = h5readatt(filename,'/data/input/model','beta');
DATA.inputs.TAU_E = h5readatt(filename,'/data/input/model','tau_e');
DATA.inputs.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
DATA.inputs.K_cB = h5readatt(filename,'/data/input/model','k_cB');
DATA.inputs.K_gB = h5readatt(filename,'/data/input/model','k_gB');
DATA.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
DATA.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
DATA.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y';
DATA.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y';
DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y';
DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y';
% Species dependent parameters
DATA.sigma = zeros(1,DATA.Na);
DATA.tau = zeros(1,DATA.Na);
DATA.q = zeros(1,DATA.Na);
DATA.K_N = zeros(1,DATA.Na);
DATA.K_T = zeros(1,DATA.Na);
spnames = {'ions','electrons'};
for ia=1:DATA.Na
spdata = ['/data/input/',spnames{ia}];
DATA.sigma(ia) = h5readatt(filename,spdata,'sigma');
DATA.tau(ia) = h5readatt(filename,spdata,'tau');
DATA.q(ia) = h5readatt(filename,spdata,'q');
DATA.K_N(ia) = h5readatt(filename,spdata,'k_N');
DATA.K_T(ia) = h5readatt(filename,spdata,'k_T');
end
DATA.spnames = spnames{1:DATA.Na};
DATA.CONAME = DATA.CO;
% Species dependent parameters
DATA.inputs.sigma = zeros(1,DATA.inputs.Na);
DATA.inputs.tau = zeros(1,DATA.inputs.Na);
DATA.inputs.q = zeros(1,DATA.inputs.Na);
DATA.inputs.K_N = zeros(1,DATA.inputs.Na);
DATA.inputs.K_T = zeros(1,DATA.inputs.Na);
spnames = {'ions','electrons'};
for ia=1:DATA.inputs.Na
spdata = ['/data/input/',spnames{ia}];
DATA.inputs.sigma(ia) = h5readatt(filename,spdata,'sigma');
DATA.inputs.tau(ia) = h5readatt(filename,spdata,'tau');
DATA.inputs.q(ia) = h5readatt(filename,spdata,'q');
DATA.inputs.K_N(ia) = h5readatt(filename,spdata,'k_N');
DATA.inputs.K_T(ia) = h5readatt(filename,spdata,'k_T');
end
DATA.inputs.spnames = spnames{1:DATA.inputs.Na};
DATA.inputs.CONAME = DATA.inputs.CO;
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 (DATA.inputs.CLOS == 0); DATA.CLOSNAME = 'Trunc.';
elseif(DATA.inputs.CLOS == 1); DATA.CLOSNAME = 'Clos. 1';
elseif(DATA.inputs.CLOS == 2); DATA.CLOSNAME = 'Clos. 2';
end
degngrad = ['P_',num2str(DATA.PMAX),'_J_',num2str(DATA.JMAX)];
degngrad = ['P_',num2str(DATA.inputs.PMAX),'_J_',num2str(DATA.inputs.JMAX)];
degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MUx]);
% if ~DATA.LINEARITY; 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
degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
DATA.inputs.CONAME,'_CLOS_',num2str(DATA.inputs.CLOS),'_mu_%0.0e'];
degngrad = sprintf(degngrad,[DATA.inputs.K_N,DATA.inputs.NU,DATA.inputs.MUx]);
% if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
resolution = [num2str(DATA.inputs.Nx),'x',num2str(DATA.inputs.Ny),'_'];
gridname = ['L_',num2str(DATA.inputs.L),'_'];
DATA.params_string = [resolution,gridname,degngrad];
end
\ No newline at end of file
matlab/plot/plot_radial_transport_and_spacetime.m
+ 8
− 67
View file @ ca8720fd
......@@ -10,12 +10,6 @@ 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(['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]);
% Ns3D = numel(DATA.Ts3D);
% [KX, KY] = meshgrid(DATA.kx, DATA.ky);
%% error estimation
DT_ = (tend-tstart)/OPTIONS.NCUT;
Qx_ee = zeros(1,OPTIONS.NCUT);
......@@ -27,41 +21,16 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee);
disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
% %% computations
%
% % Compute zonal and non zonal energies
% E_Zmode_SK = zeros(1,Ns3D);
% E_NZmode_SK = zeros(1,Ns3D);
% for it = 1:numel(DATA.Ts3D)
% E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(squeeze(f_avg_z(ikzf,1,it))).^2);
% E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0)))));
% end
% % Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP
% % 1/2 sum_p sum_j Napj^2(k=0) (avg z)
% Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj))));
% ff = trapz(DATA.z,Nipjz,5);
% E_TE = 0.5*squeeze(ff);
% % Compute electrostatic energy
% E_ES = zeros(size(DATA.Ts5D));
% bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel
% for it5D = 1:numel(DATA.Ts5D)
% [~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D)));
% for in = 1:DATA.Jmaxi
% Knphi = kernel(in-1,bi).*squeeze(trapz(DATA.z,DATA.PHI(:,:,:,it3D),3));
% Ni0n_z= squeeze(trapz(DATA.z,DATA.Nipj(1,in,:,:,:,it5D),5));
% E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z)));
% end
% end
%% Figure
clr_ = lines(20);
mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position', [500, 1000, 1000, 600])
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.params_string]; %set(gcf, 'Position', [500, 1000, 1000, 600])
FIGURE.ax1 = subplot(2,1,1,'parent',FIGURE.fig);
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
'color',clr_((DATA.Pmaxi-1)/2,:),...
'color',clr_((DATA.grids.Np-1)/2,:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
'color',clr_((DATA.Pmaxi-1)/2,:),...
'color',clr_((DATA.grids.Np-1)/2,:),...
'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
ylabel('Transport')
if(~isnan(Qx_infty_avg))
......@@ -79,24 +48,13 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
% %% Free energy
% FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig);
% yyaxis left
% plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on;
% ylabel('Entropy');%('$\epsilon_f$')
% yyaxis right
% plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$');
% ylabel('ES energy');%('$\epsilon_\phi$')
% xlim([DATA.Ts5D(1), DATA.Ts5D(end)]);
% xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
plt = @(x) mean(x(:,:,:),1);
kycut = max(DATA.ky);
kxcut = max(DATA.kx);
kycut = max(DATA.grids.ky);
kxcut = max(DATA.grids.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
OPTIONS.NAME = OPTIONS.ST_FIELD;
......@@ -104,12 +62,12 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
OPTIONS.COMP = 'avg';
OPTIONS.TIME = DATA.Ts3D;
OPTIONS.POLARPLOT = 0;
toplot = process_field(DATA,OPTIONS);
toplot = process_field_v2(DATA,OPTIONS);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
FIGURE.ax3 = subplot(2,1,2,'parent',FIGURE.fig);
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
[TY,TX] = meshgrid(DATA.grids.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
......@@ -124,21 +82,4 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
%% Zonal vs NZonal energies
% subplot(312)
% it0 = 1; itend = Ns3D;
% trange = toplot.FRAMES;
% plt1 = @(x) x;%-x(1);
% plt2 = @(x) x./max((x(:)));
% toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
% % plty = @(x) x(500:end)./max(squeeze(x(500:end)));
% yyaxis left
% % plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
% plot(plt1(DATA.Ts3D(trange)),plt2(toplot(:)),'DisplayName','Sum $A^2$');
% ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
% yyaxis right
% plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
% xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
% ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
% xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
\ No newline at end of file
matlab/plot/show_moments_spectrum.m
+ 42
− 166
View file @ ca8720fd
function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
species_name = {'i','e'}; % hard coded because this list is not output yet
Pi = DATA.Pi;
Ji = DATA.Ji;
if ~(isempty(DATA.Nipjz))
Time_ = DATA.Ts3D;
Nipj = sum(abs(DATA.Nipjz),3);
else
Time_ = DATA.Ts5D;
% Nipj = sum(sum(sum(abs(DATA.Nipj).^2,3),4),5);
Nipj = sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj,3),4),5);
end
Nipj = squeeze(Nipj);
if DATA.KIN_E
Pe = DATA.Pe;
Je = DATA.Je;
if ~(isempty(DATA.Nepjz))
Nepj = sum(abs(DATA.Nepjz),3);
else
Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
end
Nepj = squeeze(Nepj);
end
FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS];
Pa = DATA.grids.Parray;
Ja = DATA.grids.Jarray;
Time_ = DATA.Ts3D;
FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.params_string];
set(gcf, 'Position', [100 50 1000 400])
if ~OPTIONS.ST
t0 = OPTIONS.TIME(1);
t1 = OPTIONS.TIME(end);
[~,it0] = min(abs(t0-Time_));
[~,it1] = min(abs(t1-Time_));
Nipj = mean(Nipj(:,:,it0:it1),3);
if DATA.KIN_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);
ymini = min(Nipj); ymaxi = max(Nipj);
if DATA.KIN_E
Nepj = squeeze(Nepj);
ymine = min(Nepj); ymaxe = max(Nepj);
ymax = max([ymaxi ymaxe]);
ymin = min([ymini ymine]);
end
if DATA.KIN_E
subplot(121)
end
if ~OPTIONS.P2J
for ij = 1:numel(Ji)
name = ['$j=',num2str(Ji(ij)),'$'];
semilogy(Pi,Nipj(:,ij),'o-','DisplayName',name); hold on;
end
xlabel('$p$');
else
for ij = 1:numel(Ji)
name = ['$j=',num2str(Ji(ij)),'$'];
semilogy(Pi+2*Ji(ij),Nipj(:,ij),'o-','DisplayName',name); hold on;
end
xlabel('$p+2j$')
end
ylabel(['$\sum_{kx,ky}|N_i^{pj}|$']);
legend('show');
title([TITLE,' He-La ion spectrum']);
if DATA.KIN_E
subplot(122)
if ~OPTIONS.P2J
for ij = 1:numel(Je)
name = ['$j=',num2str(Je(ij)),'$'];
semilogy(Pe,Nepj(:,ij),'o-','DisplayName',name); hold on;
end
xlabel('p');
else
for ij = 1:numel(Je)
name = ['$j=',num2str(Je(ij)),'$'];
semilogy(Pe+2*Je(ij),Nepj(:,ij),'o-','DisplayName',name); hold on;
end
xlabel('$p+2j$')
end
ylabel(['$\sum_{kx,ky}|N_e^{pj}|$']);
legend('show');
title([TITLE,' He-La elec. spectrum']);
end
else
for ia = 1:DATA.inputs.Na
Napjz = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
subplot(double(DATA.inputs.Na),1,double(ia))
if OPTIONS.P2J
plotname = '$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$';
[JJ,PP] = meshgrid(Ji,Pi);
P2Ji = PP + 2*JJ;
plotname = ['$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
[JJ,PP] = meshgrid(Ja,Pa);
P2J = PP + 2*JJ;
% form an axis of velocity ordered moments
y_ = unique(P2Ji); yname = '$p+2j$';
p2j = unique(P2J); yname = '$p+2j$';
% weights to average
weights = zeros(size(y_));
weights = zeros(size(p2j));
% space time of moments amplitude wrt to p+2j degree
Ni_ST = zeros(numel(y_),numel(Time_));
Na_ST = zeros(numel(p2j),numel(Time_));
% fill the st diagramm by averaging every p+2j=n moments together
for ip = 1:numel(Pi)
for ij = 1:numel(Ji)
[~,ip2j] = min(abs(y_-(Pi(ip)+2*Ji(ij))));
Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
for ip = 1:numel(Pa)
for ij = 1:numel(Ja)
[~,ip2j] = min(abs(p2j-(Pa(ip)+2*Ja(ij))));
Na_ST(ip2j,:) = Na_ST(ip2j,:) + transpose(squeeze(Napjz(ip,ij,:)));
weights(ip2j) = weights(ip2j) + 1;
end
end
% doing the average
for ip2j = 1:numel(y_)
Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
end
if DATA.KIN_E
% same for electrons!!
[JJ,PP] = meshgrid(Je,Pe);
P2Je = PP + 2*JJ;
% form an axis of velocity ordered moments
p2je = unique(P2Je);
% weights to average
weights = zeros(size(p2je));
% space time of moments amplitude wrt to p+2j degree
Ne_ST = zeros(numel(p2je),numel(Time_));
% fill the st diagramm by averaging every p+2j=n moments together
for ip = 1:numel(Pe)
for ij = 1:numel(Je)
[~,ip2j] = min(abs(y_-(Pe(ip)+2*Je(ij))));
Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
weights(ip2j) = weights(ip2j) + 1;
end
end
% doing the average
for ip2j = 1:numel(y_)
Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
end
for ip2j = 1:numel(p2j)
Na_ST(ip2j,:) = Na_ST(ip2j,:)/weights(ip2j);
end
ticks_labels = {};
else % We just order our moments w.r.t. to the convention ..
% (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
plotname = '$\langle\sum_k |N_i^{pj}|^2\rangle_z$';
Nmoments = numel(Nipj(:,:,1)); % total number of moments
Nmoments = numel(Napjz(:,:,1)); % total number of moments
HL_deg = zeros(Nmoments,2); % list of degrees, first column Hermite, second Laguerre
im = 2;
deg = 1; % the zero degree is always here first place
......@@ -149,42 +46,34 @@ else
FOUND = 1;
while(FOUND) % As we find a pair matching the degree we retry
FOUND = 0;
for ij = 1:DATA.Jmaxi
for ip = 1:DATA.Pmaxi
if((ip-1)+2*(ij-1) == deg)
% Check if the pair is already added
check_ = HL_deg == [DATA.Pi(ip) DATA.Pi(ij)];
check_ = sum(check_(:,1) .* check_(:,2));
if ~check_ % if not add it
HL_deg(im,1) = DATA.Pi(ip);
HL_deg(im,2) = DATA.Pi(ij);
ticks_labels{im} = ['(',num2str(DATA.Pi(ip)),',',num2str(DATA.Ji(ij)),')'];
im = im + 1; FOUND = 1;
for ij = 1:DATA.grids.Nj
for ip = 1:DATA.grids.Np
if((ip-1)+2*(ij-1) == deg)
% Check if the pair is already added
check_ = HL_deg == [DATA.grids.Parray(ip) DATA.grids.Jarray(ij)];
check_ = sum(check_(:,1) .* check_(:,2));
if ~check_ % if not add it
HL_deg(im,1) = DATA.grids.Parray(ip);
HL_deg(im,2) = DATA.grids.Jarray(ij);
ticks_labels{im} = ['(',num2str(DATA.grids.Parray(ip)),',',num2str(DATA.grids.Jarray(ij)),')'];
im = im + 1; FOUND = 1;
end
end
end
end
end
end
end
% No pair found anymore, increase the degree
deg = deg + 1;
end
% form an axis of velocity ordered moments
y_ = 1:Nmoments; yname = '$(P,J)$';
% space time of moments amplitude wrt to p+2j degree
Ni_ST = zeros(Nmoments,numel(Time_));
for i_ = 1:Nmoments
Ni_ST(i_,:) = Nipj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
end
if DATA.KIN_E
p2j = 1:Nmoments; yname = '$(P,J)$';
% space time of moments amplitude wrt to p+2j degree
Ne_ST = zeros(Nmoments,numel(Time_));
Na_ST = zeros(Nmoments,numel(Time_));
for i_ = 1:Nmoments
Ne_ST = Nepj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
Na_ST(i_,:) = Napjz((HL_deg(i_,1)/DATA.grids.deltap)+1,HL_deg(i_,2)+1,:);
end
end
end
end
% plots
% scaling
if OPTIONS.NORMALIZED
......@@ -192,29 +81,16 @@ else
else
plt = @(x,i) x;
end
if DATA.KIN_E
subplot(2,1,1)
end
imagesc(Time_,y_,plt(Ni_ST,1:numel(y_)));
imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
set(gca,'YDir','normal')
xlabel('$t$'); ylabel(yname);
if ~isempty(ticks_labels)
yticks(y_);
yticks(p2j);
yticklabels(ticks_labels)
end
title(plotname)
if DATA.KIN_E
subplot(2,1,2)
imagesc(Time_,p2je,plt(Ne_ST,1:numel(y_)));
set(gca,'YDir','normal')
xlabel('$t$'); ylabel(yname)
title(plotname)
suptitle(DATA.param_title);
end
end
end
matlab/process_field_v2.m
+ 14
− 15
View file @ ca8720fd
......@@ -9,35 +9,35 @@ for i = 1:numel(OPTIONS.TIME)
end
FRAMES = unique(FRAMES);
%% Setup the plot geometry
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
directions = {'y','x','z'};
Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
[X,Y] = meshgrid(DATA.x,DATA.y);
[X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.x);
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.y);
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
[X,Y] = meshgrid(DATA.kx,DATA.ky);
[X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.kx);
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'sx'
XNAME = '$v_\parallel$'; YNAME = '$\mu$';
......@@ -112,24 +112,24 @@ switch OPTIONS.COMP
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
COMPNAME = sprintf(['y=','%2.1f'],DATA.x(i));
COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i));
else
COMPNAME = sprintf(['k_y=','%2.1f'],DATA.kx(i));
COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
COMPNAME = sprintf(['x=','%2.1f'],DATA.y(i));
COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i));
else
COMPNAME = sprintf(['k_x=','%2.1f'],DATA.ky(i));
COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
COMPNAME = sprintf(['z=','%2.1f'],DATA.grids.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
......@@ -161,7 +161,6 @@ end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
% FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
......@@ -296,7 +295,7 @@ else
tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
for iz = 1:numel(DATA.grids.z)
tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment