diff --git a/.gitignore b/.gitignore
index 622d3d8271b884ca852881cbdd72fc50677e6c18..69461096cd7858edd41a8780b7a8406fdeb00469 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,6 +23,7 @@ srcinfo
srcinfo.h
logs/
results/
+results_old/
mod/
obj/
bin/
diff --git a/matlab/LinearFit_s.m b/matlab/LinearFit_s.m
index c5348698b01f43eb8bdb35a1420f71671a6305ed..0d683f4f027e5a34c9232e92ecf21508153b467f 100644
--- a/matlab/LinearFit_s.m
+++ b/matlab/LinearFit_s.m
@@ -22,7 +22,7 @@ function [gamma,fit] = LinearFit_s(time,Na00abs, tstart, tend)
Na00absshifted = circshift(Na00abs,-1); % ... shift by -1 the time position
- gammaoft = log(Na00absshifted(1:end-1)./Na00abs(1:end-1))./transpose(diff(time)); % ... evaluate growth rate
+ gammaoft = log(Na00absshifted(1:end-1)./Na00abs(1:end-1))./squeeze(diff(time)); % ... evaluate growth rate
% Get gamma
gamma = mean(gammaoft(Tstart_ind:Tend_ind-1)); % ... take the mean of gamma over the time window
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index f4bbf44cbd9471fe89c0ef7449fd4c9698c45bb8..826ea6def012776a647ec57e97ccd545e6f2ed0b 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -1,9 +1,5 @@
title1 = GIFNAME;
-FIGDIR = ['../results/',BASIC.SIMID,'/'];
-if ~exist(FIGDIR, 'dir')
- mkdir(FIGDIR)
-end
-
+FIGDIR = BASIC.RESDIR;
GIFNAME = [FIGDIR, GIFNAME,'.gif'];
% Set colormap boundaries
diff --git a/matlab/create_gif_1D.m b/matlab/create_gif_1D.m
index f0a0e8e8f2b09b78c661faf3e422615bf2479462..5f7369fafabe42a1fca89c2384a520d404db360c 100644
--- a/matlab/create_gif_1D.m
+++ b/matlab/create_gif_1D.m
@@ -1,5 +1,5 @@
title1 = GIFNAME;
-FIGDIR = ['../results/',BASIC.SIMID,'/'];
+FIGDIR = BASIC.RESDIR;
if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
diff --git a/matlab/save_figure.m b/matlab/save_figure.m
index 9fae159ba479532926bd5f2df9f5bb50e17f1191..e71ee646951dc1664ee43b5f009d07f451298b36 100644
--- a/matlab/save_figure.m
+++ b/matlab/save_figure.m
@@ -1,10 +1,5 @@
%% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME)
% and parameters
-FIGDIR = ['../results/', BASIC.SIMID,'/'];
-if ~exist(FIGDIR, 'dir')
- mkdir(FIGDIR)
-end
-
-FIGNAME = [FIGDIR, FIGNAME,FMT];
+FIGNAME = [BASIC.RESDIR, FIGNAME,FMT];
saveas(fig,FIGNAME);
disp(['Figure saved @ : ',FIGNAME])
\ No newline at end of file
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 8d9fd05acafc820db7d2f90640cfd5fae90411f2..f9f3fa748ffe478a1f00dc77b41102fbc17f0638 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -1,6 +1,6 @@
function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC)
% Write the input script "fort.90" with desired parameters
-INPUT = '../wk/fort.90';
+INPUT = 'fort.90';
fid = fopen(INPUT,'wt');
fprintf(fid,'&BASIC\n');
@@ -73,4 +73,5 @@ fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
fprintf(fid,'/');
fclose(fid);
+system(['cp fort.90 ',BASIC.RESDIR,'/.']);
end
diff --git a/wk/analysis_1D.m b/wk/analysis_1D.m
index 31d07cd51baaf996a969cde9632419637d5b7e2f..9076e8565eaacbc69116971b7f28bf3825f0267b 100644
--- a/wk/analysis_1D.m
+++ b/wk/analysis_1D.m
@@ -1,16 +1,13 @@
default_plots_options
-%% load results
+%% load resulTs2D
JOBNUM = 00;
-filename = [BASIC.SIMID,'_','%.2d.h5'];
-filename = sprintf(filename,JOBNUM); disp(['Analysing ',filename])
-[Nipj, p_, j_, kr, kz, Ts, dt] = load_5D_data(filename, 'moments_i');
-Nepj = load_5D_data(filename, 'moments_e');
-Ni00 = squeeze(Nipj(1,1,:,:,:));
-Ne00 = squeeze(Nepj(1,1,:,:,:));
-PH = squeeze(load_2D_data(filename, 'phi'));
-Ts = Ts';
-Ns = numel(Ts);
-dt_samp = mean(diff(Ts));
+load_results
+Ni00 = squeeze(Ni00);
+Ne00 = squeeze(Ne00);
+PHI = squeeze(PHI);
+Ts2D = Ts2D';
+Ns = numel(Ts2D);
+dt_samp = mean(diff(Ts2D));
% Build grids
Nkr = numel(kr); Nkz = numel(kz);
[KZ,KR] = meshgrid(kz,kr);
@@ -29,8 +26,8 @@ ne00 = zeros(Nz, Ns);
ni00 = zeros(Nz, Ns);
phi = zeros(Nz, Ns);
-for it = 1:numel(PH(1,:))
- NE_ = Ne00(:,it); NN_ = Ni00(:,it); PH_ = PH(:,it);
+for it = 1:numel(Ts2D)
+ NE_ = Ne00(:,it); NN_ = Ni00(:,it); PH_ = PHI(:,it);
F_ = (ifft((NE_),Nz));
ne00(:,it)= real(fftshift(F_));
F_ = (ifft((NN_),Nz));
@@ -41,9 +38,7 @@ end
% Post processing
ne_00 = zeros(1,Ns); % Time evolution of ne(r,z) at origin
-Ne_gm = zeros(1,Ns); % Time evolution of Ne(k) max gamma (max real)
ni_00 = zeros(1,Ns); % .
-Ni_gm = zeros(1,Ns); % .
phi_00 = zeros(1,Ns); % .
E_pot = zeros(1,Ns); % Potential energy n^2
E_kin = zeros(1,Ns); % Kinetic energy grad(phi)^2
@@ -52,14 +47,12 @@ CFL = zeros(1,Ns); % CFL time step
Ddz = 1i*kz;
[~,iz0] = min(abs(z)); % index of z==0
[~,ikz1] = min(abs(kz-round(1/dkz)*dkz)); % index of kz==1
-for it = 1:numel(PH(1,:))
- NE_ = squeeze(Ne00(:,it)); NN_ = squeeze(Ni00(:,it)); PH_ = squeeze(PH(:,it));
+for it = 1:numel(Ts2D)
+ NE_ = squeeze(Ne00(:,it)); NN_ = squeeze(Ni00(:,it)); PH_ = squeeze(PHI(:,it));
ne_00(it) = ne00(iz0,it);
- Ne_gm(it) = max(real(Ne00(:,it)));
ni_00(it) = ni00(iz0,it);
- Ni_gm(it) = max(real(Ni00(:,it)));
phi_00(it) = phi(iz0,it);
@@ -67,79 +60,86 @@ for it = 1:numel(PH(1,:))
E_kin(it) = pi/Lz*sum(abs(Ddz.*PH_).^2)/Nkz;
ExB(it) = max(abs(phi(3:end,it)-phi(1:end-2,it))'/(2*dz));
end
-E_kin_KZ = abs(Ddz.*PH(:,it)).^2;
-dEdt = diff(E_pot+E_kin)./diff(Ts);
+E_kin_KZ = abs(Ddz.*PHI(:,it)).^2;
+dEdt = diff(E_pot+E_kin)./diff(Ts2D);
-% Growth rate
+%% Compute growth rate
gamma_Ne = zeros(1,Nkz);
gamma_Ni = zeros(1,Nkz);
gamma_PH = zeros(1,Nkz);
+tend = Ts2D(end); tstart = 0.6*tend;
+plt = @(x) squeeze(abs(x));
for ikz = 1:Nkz
- gamma_Ne(ikz) = real(LinearFit_s(Ts,Ne00(ikz,:),-1,-1));
- gamma_Ni(ikz) = real(LinearFit_s(Ts,Ni00(ikz,:),-1,-1));
- gamma_PH(ikz) = real(LinearFit_s(Ts,PH(ikz,:),-1,-1));
+ gamma_Ne(ikz) = LinearFit_s(Ts2D,plt(Ne00(ikz,:)),tstart,tend);
+ gamma_Ni(ikz) = LinearFit_s(Ts2D,plt(Ni00(ikz,:)),tstart,tend);
+ gamma_PH(ikz) = LinearFit_s(Ts2D,plt(PHI(ikz,:)),tstart,tend);
end
-gamma_Ne = gamma_Ne .* (gamma_Ne>=0.0);
-gamma_Ni = gamma_Ni .* (gamma_Ni>=0.0);
-gamma_PH = gamma_PH .* (gamma_PH>=0.0);
-%% PLOTS
+gamma_Ne = real(gamma_Ne .* (gamma_Ne>=0.0));
+gamma_Ni = real(gamma_Ni .* (gamma_Ni>=0.0));
+gamma_PH = real(gamma_PH .* (gamma_PH>=0.0));
+
+%% PLOTs2D
+if 0
%% Time evolutions
fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
subplot(221)
- semilogy(Ts,abs(ne_00),'-','DisplayName','$n_e^{00}$'); hold on;
- semilogy(Ts,abs(ni_00),'-','DisplayName','$n_i^{00}$');
+ semilogy(Ts2D,abs(ne_00),'-','DisplayName','$n_e^{00}$'); hold on;
+ semilogy(Ts2D,abs(ni_00),'-','DisplayName','$n_i^{00}$');
grid on; xlabel('$t$'); ylabel('$|n_a(x=0,y=0)|$');
subplot(222)
- semilogy(Ts,abs(Ni_gm),'-','DisplayName','$\phi$')
+ semilogy(Ts2D,abs(Ni_gm),'-','DisplayName','$\phi$')
grid on; xlabel('$t$'); ylabel('$|\tilde n(k_r\approx 1,k_z\approx 1)|$');
subplot(223)
- semilogy(Ts,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
+ semilogy(Ts2D,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
hold on;
grid on; xlabel('$t$'); ylabel('$E$'); legend('show');
FMT = '.fig'; save_figure
-
+end
+if 1
%% Growth rate
fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
+plt = @(x) circshift(x,Nkz/2-1);
subplot(221)
- plot(kz,gamma_Ne,'-'); hold on;
+ plot(plt(kz),plt(gamma_Ne),'-'); hold on; xlim([0.0,max(kz)+dkz]);
grid on; xlabel('$k_z$'); ylabel('$\gamma(N_e^{00})$');
subplot(222)
- plot(kz,gamma_Ni,'-'); hold on;
+ plot(plt(kz),plt(gamma_Ni),'-'); hold on; xlim([0.0,max(kz)+dkz]);
grid on; xlabel('$k_z$'); ylabel('$\gamma(N_i^{00})$');
subplot(223)
- plot(kz,gamma_PH,'-'); hold on;
+ plot(plt(kz),plt(gamma_PH),'-'); hold on; xlim([0.0,max(kz)+dkz]);
grid on; xlabel('$k_z$'); ylabel('$\gamma(\tilde\phi)$'); legend('show');
FMT = '.fig'; save_figure
-
+end
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0 = 0;
skip_ = 1;
DELAY = 0.01*skip_;
-FRAMES = floor(t0/dt_samp)+1:skip_:numel(Ts);
+FRAMES = floor(t0/dt_samp)+1:skip_:numel(Ts2D);
+linestyle = 'o-.';
if 0
%% Density electron
GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)];
-FIELD = real(ne00); X = z; T = Ts;
+FIELD = real(ne00); X = z; T = Ts2D;
FIELDNAME = '$n_e^{00}/\max(n_e^{00})$'; XNAME = '$z$';
XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
create_gif_1D
%% Density ion
GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)];
-FIELD = real(ni00); X = z; T = Ts;
+FIELD = real(ni00); X = z; T = Ts2D;
FIELDNAME = '$n_i^{00}/\max(n_i^{00})$'; XNAME = '$z$';
XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
create_gif_1D
%% Phi
GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];
-FIELD = real(phi); X = z; T = Ts;
+FIELD = real(phi); X = z; T = Ts2D;
FIELDNAME = '$\phi/\max(\phi)$'; XNAME = '$z$';
XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
create_gif_1D
%% Density electron frequency
GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)];
-FIELD = (abs(Ni00)); X = (kz); T = Ts;
+FIELD = (abs(Ni00)); X = (kz); T = Ts2D;
FIELDNAME = '$|N_i^{00}|$'; XNAME = '$k_z$';
-XMIN = -0.5; XMAX = max(kz)+.5; YMIN = -0.1; YMAX = 1.1;
+XMIN = min(kz)-0.5; XMAX = max(kz)+.5; YMIN = -0.1; YMAX = 1.1;
create_gif_1D
end
@@ -147,7 +147,7 @@ if 0
%% Space-Time diagrams at r = 0
plt = @(x) real(x);
fig = figure; FIGNAME = ['z_space_time_diag',sprintf('_%.2d',JOBNUM)];
- [TY,TX] = meshgrid(Ts,z);
+ [TY,TX] = meshgrid(Ts2D,z);
subplot(221)% density
pclr = pcolor(TX,TY,(plt(ne00))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_e^{00}$');
@@ -162,7 +162,7 @@ FMT = '.fig'; save_figure
%% Space-Time diagrams at kr=0
plt = @(x) abs(x);
fig = figure; FIGNAME = ['kz_space_time_diag',sprintf('_%.2d',JOBNUM)];
- [TY,TX] = meshgrid(Ts,kz);
+ [TY,TX] = meshgrid(Ts2D,kz);
subplot(221)% density
pclr = pcolor(TX,TY,(plt(Ne00))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$kz$'); ylabel('$t$'); title('$\textrm{Re}(N_e^{00})|_{k_r=0}$');
@@ -183,41 +183,51 @@ if 0
[~,ik20] = min(abs(kz-2.00));
fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
- subplot(221); plt = @(x) (abs(x));
- semilogy(Ts,plt(PH(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts,plt(PH(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts,plt(PH(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts,plt(PH(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$\hat\phi$'); legend('show');
- semilogy(Ts,abs(PH(ik05,end)).*exp(gamma_PH(ik05).*(Ts-Ts(end))),'--k')
- semilogy(Ts,abs(PH(ik10,end)).*exp(gamma_PH(ik10).*(Ts-Ts(end))),'--k')
- subplot(222); plt = @(x) (abs(x));
- semilogy(Ts,plt(Ni00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts,plt(Ni00(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts,plt(Ni00(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts,plt(Ni00(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show');
- subplot(223); plt = @(x) (abs(x));
- semilogy(Ts,plt(Ne00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts,plt(Ne00(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts,plt(Ne00(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts,plt(Ne00(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$\hat n_e^{00}$'); legend('show');
+ subplot(221); plt = @(x) abs(squeeze(x));
+ semilogy(Ts2D,plt(PHI(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts2D,plt(PHI(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts2D,plt(PHI(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts2D,plt(PHI(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$|\tilde\phi|$'); legend('show');
+ semilogy(Ts2D,plt(PHI(ik05,end)).*exp(gamma_PH(ik05).*(Ts2D-Ts2D(end))),'--k')
+ semilogy(Ts2D,plt(PHI(ik10,end)).*exp(gamma_PH(ik10).*(Ts2D-Ts2D(end))),'--k')
+ plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
+ plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
+ subplot(222); plt = @(x) squeeze(abs(x));
+ semilogy(Ts2D,plt(Ni00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts2D,plt(Ni00(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts2D,plt(Ni00(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts2D,plt(Ni00(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$|N_i^{00}|$'); legend('show');
+ semilogy(Ts2D,plt(Ni00(ik05,end)).*exp(gamma_Ni(ik05).*(Ts2D-Ts2D(end))),'--k')
+ semilogy(Ts2D,plt(Ni00(ik10,end)).*exp(gamma_Ni(ik10).*(Ts2D-Ts2D(end))),'--k')
+ plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
+ plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
+ subplot(223); plt = @(x) squeeze(abs(x));
+ semilogy(Ts2D,plt(Ne00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts2D,plt(Ne00(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts2D,plt(Ne00(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts2D,plt(Ne00(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$|N_e^{00}|$'); legend('show');
+ semilogy(Ts2D,plt(Ne00(ik05,end)).*exp(gamma_Ne(ik05).*(Ts2D-Ts2D(end))),'--k')
+ semilogy(Ts2D,plt(Ne00(ik10,end)).*exp(gamma_Ne(ik10).*(Ts2D-Ts2D(end))),'--k')
+ plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
+ plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
FMT = '.fig'; save_figure
end
if 0
%% Show frame
-it = min(70,numel(Ts));
+it = min(70,numel(Ts2D));
fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
subplot(221); plt = @(x) (abs(x));
plot(kz,plt(PH(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat\phi$');
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat\phi$');
subplot(222); plt = @(x) (abs(x));
plot(kz,plt(Ni00(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat n_i^{00}$');
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat n_i^{00}$');
subplot(223); plt = @(x) (abs(x));
plot(kz,plt(Ne00(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat n_e^{00}$');
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat n_e^{00}$');
FMT = '.fig'; save_figure
end
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 5c3fd115d7e0b35d87e8f2885397c52dce34934e..1ffdb66d1705c94f5184869cc8f5dcfb5d9f0693 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,5 +1,5 @@
%% Load results
-load_results
+compile_results
%% Retrieving max polynomial degree and sampling info
Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -76,10 +76,10 @@ for it = 1:numel(Ts2D) % Loop over 2D arrays
E_pot(it) = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id
E_kin(it) = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr;
ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
- Flux_ri(it) = sum(sum(ni00(:,:,it).*drphi(:,:,it)))*dr*dz;
- Flux_zi(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz;
- Flux_re(it) = sum(sum(ne00(:,:,it).*drphi(:,:,it)))*dr*dz;
- Flux_ze(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz;
+ Flux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
+ Flux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
+ Flux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
+ Flux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
end
E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
@@ -88,18 +88,18 @@ dEdt = diff(E_pot+E_kin)./dt2D;
for it = 1:numel(Ts5D) % Loop over 5D arrays
NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
- Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4);
- Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4);
+ Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
+ Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
if strcmp(OUTPUTS.write_non_lin,'.true.')
- Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4);
- Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))));
- Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4);
- Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))));
+ Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
+ Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz;
+ Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
+ Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz;
end
end
-%% Growth rate
+%% Compute growth rate
disp('- growth rate')
tend = Ts2D(end); tstart = 0.6*tend;
g_ = zeros(Nkr,Nkz);
@@ -114,6 +114,7 @@ gkr0kz_Ni00 = real(g_(ikr0KH,:));
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Plots')
+FMT = '.fig';
%% Time evolutions
fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
subplot(221);
@@ -151,7 +152,7 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
end
grid on; xlabel('$t$'); ylabel('$S$'); %legend('show');
else
- %% Growth rate
+%% Growth rate
subplot(224)
[~,ikr0KH] = min(abs(kr-KR0KH));
plot(kz(1:ikr0KH)/kr(ikr0KH),gkr0kz_Ni00(1:ikr0KH)/(KR0KH*A0KH),...
@@ -161,33 +162,63 @@ else
end
save_figure
+%%
+if 1
+%% Show frame in real space
+tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
+fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)];
+ subplot(221); plt = @(x) (((x)));
+ pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_e^{00}|$');
+ subplot(222); plt = @(x) ((x));
+ pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_i^{00}|$');
+ subplot(223); plt = @(x) ((x));
+ pclr = pcolor((RR),(ZZ),plt(phi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat\phi|$');
+if strcmp(OUTPUTS.write_non_lin,'.true.')
+ subplot(224); plt = @(x) fftshift((abs(x)),2);
+ pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it5D))); legend('$|S_i^{00}|$');
+end
+save_figure
+end
+
%%
t0 = 0;
skip_ = 1;
DELAY = 0.01*skip_;
FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
-if 0
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if 0
%% Density ion
GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D;
FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
create_gif
+end
+if 0
%% Density electron
GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D;
FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$';
create_gif
+end
+if 0
%% Phi
GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D;
FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$';
create_gif
+end
+if 0
%% Density ion frequency
GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
create_gif
+end
+if 0
%% Density ion frequency @ kr = 0
GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.';
@@ -290,28 +321,6 @@ fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)];
save_figure
end
-%%
-if 0
-%% Show frame in real space
-tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
-fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)];
- subplot(221); plt = @(x) (((x)));
- pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_e^{00}|$');
- subplot(222); plt = @(x) ((x));
- pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_i^{00}|$');
- subplot(223); plt = @(x) ((x));
- pclr = pcolor((RR),(ZZ),plt(phi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat\phi|$');
-if strcmp(OUTPUTS.write_non_lin,'.true.')
- subplot(224); plt = @(x) fftshift((abs(x)),2);
- pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it5D))); legend('$|S_i^{00}|$');
-end
-save_figure
-end
-
%%
if 0
%% Show frame in kspace
@@ -333,12 +342,3 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
end
save_figure
end
-
-
-if 0
-%% Check time evolution of higher moments
-p = 0; j = 0;
-if strcmp(OUTPUTS.write_moments,'.true.')
-
-end
-end
diff --git a/wk/load_results.m b/wk/load_results.m
index a8c70c175850c74dd78784c265cbfdba5de289a1..5fccf08ef6a800454e2452f411f99ca39804dca3 100644
--- a/wk/load_results.m
+++ b/wk/load_results.m
@@ -1,8 +1,6 @@
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-JOBNUM = 00;
-FMT = '.fig';
-filename = [BASIC.SIMID,'_','%.2d.h5'];
-filename = sprintf(filename,JOBNUM); disp(['Loading ',filename])
+filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
+disp(['Loading ',filename])
% Loading from output file
% CPUTIME = h5readatt(filename,'/data/input','cpu_time');
if strcmp(OUTPUTS.write_moments,'.true.')
diff --git a/wk/parameters_ZP_linear.m b/wk/parameters_ZP_linear.m
deleted file mode 100644
index ceebcb0cd6a7e03b2c9b905c436de2db9e28c988..0000000000000000000000000000000000000000
--- a/wk/parameters_ZP_linear.m
+++ /dev/null
@@ -1,44 +0,0 @@
-clear all;
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 1e-2; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-ETAB = 0.5; % Magnetic gradient
-ETAN = 1.0; % Density gradient
-ETAT = 0.0; % Temperature gradient
-MU = 1e-4; % Hyper diffusivity coefficient
-LAMBDAD = 0.0;
-NOISE0 = 5.0e-5;
-%% GRID PARAMETERS
-N = 64; % Frequency gridpoints (Nkr = N/2)
-L = 100; % Size of the squared frequency domain
-KREQ0 = 1; % put kr = 0
-PMAXE = 02; % Highest electron Hermite polynomial degree
-JMAXE = 01; % Highest '' Laguerre ''
-PMAXI = 02; % Highest ion Hermite polynomial degree
-JMAXI = 01; % Highest '' Laguerre ''
-KPAR = 0.0; % Parellel wave vector component
-%% TIME PARAMETERS
-TMAX = 300; % Maximal time unit
-DT = 1e-2; % Time step
-SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1; % Sampling per time unit for 5D arrays
-RESTART = 0; % To restart from last checkpoint
-JOB2LOAD= 00;
-%% OPTIONS
-SIMID = 'ZP_'; % Name of the simulation
-NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
-CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
-KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-NO_E = 0; % Remove electrons dynamic
-% DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
-
-
-setup
diff --git a/wk/setup.m b/wk/setup.m
index d47f3b11e46afd1a60fa108b2bd471d7d4b7ffda..8b13ed6aae7396edc4800d6a229039038dd7b170 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -1,4 +1,5 @@
%% ________________________________________________________________________
+SIMDIR = ['../results/',SIMID,'/'];
% Grid parameters
GRID.pmaxe = PMAXE; % Electron Hermite moments
GRID.jmaxe = JMAXE; % Electron Laguerre moments
@@ -53,19 +54,20 @@ INITIAL.iemat_file = ...
'_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
% Naming and creating input file
-if (CO == 0); CONAME = '';
+if (CO == 0); CONAME = 'LB';
elseif(CO == -1); CONAME = 'FC';
elseif(CO == -2); CONAME = 'DG';
end
-params = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
+degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
'_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI),...
'_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_nu_%0.0e_',...
- CONAME,'_mu_%0.0e_'];
-params = sprintf(params,[NU,MU]);
-if ~NON_LIN; params = ['lin_',params]; end
-resolution = ['_',num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
+ CONAME,'_mu_%0.0e'];
+degngrad = sprintf(degngrad,[NU,MU]);
+if ~NON_LIN; degngrad = ['lin_',degngrad]; end
+resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
gridname = ['L_',num2str(L),'_'];
-BASIC.SIMID = [SIMID,resolution,gridname,params];
+PARAMS = [resolution,gridname,degngrad];
+BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
BASIC.nrun = 1e8;
BASIC.dt = DT;
BASIC.tmax = TMAX; %time normalized to 1/omega_pe
@@ -81,16 +83,24 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_non_lin = OUTPUTS.write_moments;
if NON_LIN == 0; OUTPUTS.write_non_lin = '.false.'; end;
OUTPUTS.write_doubleprecision = '.true.';
-OUTPUTS.resfile0 = ['''',BASIC.SIMID,''''];
-OUTPUTS.rstfile0 = ['''','../checkpoint/cp_',BASIC.SIMID,''''];
+OUTPUTS.resfile0 = ['''',BASIC.RESDIR,'outputs',''''];
+OUTPUTS.rstfile0 = ['''',BASIC.RESDIR,'checkpoint',''''];
OUTPUTS.job2load = JOB2LOAD;
-%% Compile and write input file
+%% Create directories
+if ~exist(SIMDIR, 'dir')
+ mkdir(SIMDIR)
+end
+if ~exist(BASIC.RESDIR, 'dir')
+mkdir(BASIC.RESDIR)
+end
+%% Compile and WRITE input file
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
MAKE = 'cd ..; make; cd wk';
system(MAKE);
%%
-disp(['Set up ', BASIC.SIMID]);
+disp(['Set up ',SIMID]);
+disp([resolution,gridname,degngrad]);
if RESTART
disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else
disp(['- starting from T = 0']);