diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index d0beda9da004f60fa15035417bb9a6db523e1c74..bb7cea95cc0fcc5aa75ac597aaf4b338b8d63fa7 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -1,4 +1,4 @@
-function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT)
+function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, OPTIONS)
% Remove AA part
if DATA.Nx > 1
@@ -11,8 +11,8 @@ ikynz = logical(ikynz .* (DATA.ky > 0));
phi = DATA.PHI(ikynz,ikxnz,:,:);
t = DATA.Ts3D;
-[~,its] = min(abs(t-TRANGE(1)));
-[~,ite] = min(abs(t-TRANGE(end)));
+[~,its] = min(abs(t-OPTIONS.TRANGE(1)));
+[~,ite] = min(abs(t-OPTIONS.TRANGE(end)));
w_ky = zeros(sum(ikynz),ite-its);
ce = zeros(sum(ikynz),ite-its);
@@ -34,7 +34,7 @@ for it = its+1:ite
end
[kys, Is] = sort(DATA.ky(ikynz));
-linear_gr.trange = t(its:ite);
+linear_gr.OPTIONS.TRANGE = t(its:ite);
linear_gr.g_ky = real(w_ky(Is,:));
linear_gr.w_ky = imag(w_ky(Is,:));
linear_gr.ce = abs(ce);
@@ -44,20 +44,29 @@ linear_gr.avg_g = mean(real(w_ky(Is,:)),2);
linear_gr.std_w = std (imag(w_ky(Is,:)),[],2);
linear_gr.avg_w = mean(imag(w_ky(Is,:)),2);
-if PLOT >0
+if OPTIONS.NPLOTS >0
figure
-if PLOT > 1
+if OPTIONS.NPLOTS > 1
subplot(1,2,1)
end
+ x_ = linear_gr.ky;
+ plt = @(x) x;
+ OVERK = '';
+ if OPTIONS.GOK == 1
+ plt = @(x) x./x_;
+ OVERK = '/$k_y$';
+ elseif OPTIONS.GOK == 2
+ plt = @(x) x.^2./x_.^3;
+ OVERK = '/$k_y$';
+ end
+ plot(x_,plt(linear_gr.g_ky(:,end)),'-o','DisplayName',['$Re(\omega_{k_y})$',OVERK]); hold on;
+ plot(x_,plt(linear_gr.w_ky(:,end)),'--*','DisplayName',['$Im(\omega_{k_y})$',OVERK]); hold on;
+ plot(x_,plt(linear_gr.ce (:,end)),'x','DisplayName',['$\epsilon$',OVERK]); hold on;
- plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
- plot(linear_gr.ky,linear_gr.w_ky(:,end),'--*','DisplayName','$Im(\omega_{k_y})$'); hold on;
- plot(linear_gr.ky,linear_gr.ce (:,end),'x','DisplayName','$\epsilon$'); hold on;
-
- errorbar(linear_gr.ky,linear_gr.avg_g,linear_gr.std_g,...
+ errorbar(linear_gr.ky,plt(linear_gr.avg_g),plt(linear_gr.std_g),...
'-o','DisplayName','$Re(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
- errorbar(linear_gr.ky,linear_gr.avg_w,linear_gr.std_w,...
+ errorbar(linear_gr.ky,plt(linear_gr.avg_w),plt(linear_gr.std_w),...
'--*','DisplayName','$Im(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
@@ -66,10 +75,10 @@ end
legend('show');
title(DATA.param_title);
-if PLOT > 1
- if PLOT == 2
+if OPTIONS.NPLOTS > 1
+ if OPTIONS.NPLOTS == 2
subplot(1,2,2)
- elseif PLOT == 3
+ elseif OPTIONS.NPLOTS == 3
subplot(2,2,2)
end
plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on;
@@ -77,7 +86,7 @@ if PLOT > 1
xlabel('t'); ylabel('$\gamma(t),\omega(t)$'); xlim([DATA.Ts3D(1) DATA.Ts3D(end)]);
end
-if PLOT > 2
+if OPTIONS.NPLOTS > 2
xlabel([]); xticks([]);
subplot(2,2,4)
[~,ikx0] = min(abs(DATA.kx));
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 1363616c9916aef89192fbca444ac0504ea2dfef..34df7362df15d32c724bb8be77521a9d3eeb3445 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,14 +1,14 @@
% Options
-SHOW_FILM = 1;
+SHOW_FILM = 0;
field2plot ='phi';
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME = 15000; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME = 1000; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
Tfin = 100;
dt_ = 0.1;
Nstep = ceil(Tfin/dt_);
% Init tracers
-Np = 20; %number of tracers
+Np = 200; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
@@ -182,7 +182,7 @@ while(t_<Tfin && it <= Nstep)
% push the particle
% q = sign(-u___(3));
- q = -u___(3);
+ q = 1;%-u___(3);
% q =1;
x_ = x_ + dt_*u___(1)*q;
y_ = y_ + dt_*u___(2)*q;
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 0a4fee17d6e512396b661d009476b1d94a54a42f..37b545481408c1d8c53295ba825db9d53acea4f3 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -8,9 +8,24 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
cstart= h5readatt(filename,'/data/input','start_iframe3d');
data = zeros(numel(ky),numel(kx),numel(z),numel(time));
- for it = 1:numel(time)
- tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
- data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ [KX,KY] = meshgrid(kx,ky);
+
+ switch variablename
+ case 'vEx'
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = 1i.*KY.*(tmp.real + 1i * tmp.imaginary);
+ end
+ case 'vEy'
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = -1i.*KX.*(tmp.real + 1i * tmp.imaginary);
+ end
+ otherwise
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
end
end
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 544f4799cf42867510c3dad5a692dd71b9823f9c..32653c85cff97c84df96cf0d0b3065e3e29a2d8d 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -195,10 +195,14 @@ resdir ='';
% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan';
% resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan';
% % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
-resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
+% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
-resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
+
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
+resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
JOBNUMMIN = 00; JOBNUMMAX = 10;
resdir = ['results/',resdir];
diff --git a/wk/lin_EPY.m b/wk/lin_EPY.m
index 429c1326ac7e1033a481c0281fdbe7f9fa620cb3..b49201d025f2e9f4872fcc826be6d1d0040bc6c0 100644
--- a/wk/lin_EPY.m
+++ b/wk/lin_EPY.m
@@ -4,17 +4,19 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_EPY'; % Name of the simulation
+SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
PROGDIR = '/home/ahoffman/gyacomo/';
-EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 2.2; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
@@ -24,14 +26,14 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
+P = 2;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 2; % real space x-gridpoints
-NY = 100; % '' y-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 120;%2*pi/0.05; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
@@ -44,7 +46,7 @@ SHEAR = 0.8; % magnetic shear
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 50; % 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
@@ -56,7 +58,7 @@ JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'LD';
+CO = 'SG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -93,7 +95,8 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
@@ -107,9 +110,10 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
%% Short analysis
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+lg = compute_fluxtube_growth_rate(data,options);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index ab514b2030232de1869342609a0e144933ea0dd2..3ebe622d264ab21cb5e9b2f822abf10c897aed18 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -4,7 +4,7 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_ITG'; % Name of the simulation
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/gyacomo/';
diff --git a/wk/save_iFFT.m b/wk/save_iFFT.m
index d429b62274e33200d342a9d95540c2ef490ec0c6..d9b3dc9c1912c1aee5918ec63cab20506d3961e9 100644
--- a/wk/save_iFFT.m
+++ b/wk/save_iFFT.m
@@ -3,7 +3,7 @@ simdir = '/misc/gyacomo_outputs/results/Zpinch_rerun';
resolu = 'UHD_512x256x2x1';
output = 'outputs_01.h5';
filename = [simdir,'/',resolu,'/',output];
-fieldname = 'Ni00';
+fieldname = 'vEy';
%% OUTPUT
outdir = '.';