From 089bc1ca8e5bed15827b22b24b26228ccdf0976b Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 6 Jul 2021 10:39:05 +0200
Subject: [PATCH] matlab scripts update
---
matlab/bluewhitered.m | 116 ++++++++++++
matlab/compile_results.m | 40 +++-
matlab/create_gif.m | 4 +-
matlab/create_mov.m | 3 +-
matlab/kernel.m | 2 +-
matlab/load_params.m | 7 +-
matlab/photomaton.m | 110 -----------
matlab/plot_kernels.m | 46 +++++
matlab/profiler.m | 28 +--
matlab/setup.m | 80 ++------
matlab/write_fort90.m | 3 -
wk/ZF_fourier_analysis.m | 24 ++-
wk/analysis_2D.m | 279 ++++++++++++++++++----------
wk/continue_multiple_runs_marconi.m | 27 ++-
wk/linear_study.m | 56 +++---
wk/load_multiple_outputs_marconi.m | 9 +-
wk/local_run.m | 22 +--
wk/marconi_run.m | 33 ++--
wk/open_figure_script.m | 20 +-
wk/photomaton.m | 90 +++++++++
20 files changed, 605 insertions(+), 394 deletions(-)
create mode 100644 matlab/bluewhitered.m
delete mode 100644 matlab/photomaton.m
create mode 100644 matlab/plot_kernels.m
create mode 100644 wk/photomaton.m
diff --git a/matlab/bluewhitered.m b/matlab/bluewhitered.m
new file mode 100644
index 00000000..b4a00b54
--- /dev/null
+++ b/matlab/bluewhitered.m
@@ -0,0 +1,116 @@
+function newmap = bluewhitered(m)
+%BLUEWHITERED Blue, white, and red color map.
+% BLUEWHITERED(M) returns an M-by-3 matrix containing a blue to white
+% to red colormap, with white corresponding to the CAXIS value closest
+% to zero. This colormap is most useful for images and surface plots
+% with positive and negative values. BLUEWHITERED, by itself, is the
+% same length as the current colormap.
+%
+% Examples:
+% ------------------------------
+% figure
+% imagesc(peaks(250));
+% colormap(bluewhitered(256)), colorbar
+%
+% figure
+% imagesc(peaks(250), [0 8])
+% colormap(bluewhitered), colorbar
+%
+% figure
+% imagesc(peaks(250), [-6 0])
+% colormap(bluewhitered), colorbar
+%
+% figure
+% surf(peaks)
+% colormap(bluewhitered)
+% axis tight
+%
+% See also HSV, HOT, COOL, BONE, COPPER, PINK, FLAG,
+% COLORMAP, RGBPLOT.
+if nargin < 1
+ m = size(get(gcf,'colormap'),1);
+end
+bottom = [0 0 0.5];
+botmiddle = [0 0.5 1];
+middle = [1 1 1];
+topmiddle = [1 0 0];
+top = [0.5 0 0];
+% Find middle
+lims = get(gca, 'CLim');
+% Find ratio of negative to positive
+if (lims(1) < 0) & (lims(2) > 0)
+ % It has both negative and positive
+ % Find ratio of negative to positive
+ ratio = abs(lims(1)) / (abs(lims(1)) + lims(2));
+ neglen = round(m*ratio);
+ poslen = m - neglen;
+
+ % Just negative
+ new = [bottom; botmiddle; middle];
+ len = length(new);
+ oldsteps = linspace(0, 1, len);
+ newsteps = linspace(0, 1, neglen);
+ newmap1 = zeros(neglen, 3);
+
+ for i=1:3
+ % Interpolate over RGB spaces of colormap
+ newmap1(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+ end
+
+ % Just positive
+ new = [middle; topmiddle; top];
+ len = length(new);
+ oldsteps = linspace(0, 1, len);
+ newsteps = linspace(0, 1, poslen);
+ newmap = zeros(poslen, 3);
+
+ for i=1:3
+ % Interpolate over RGB spaces of colormap
+ newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+ end
+
+ % And put 'em together
+ newmap = [newmap1; newmap];
+
+elseif lims(1) >= 0
+ % Just positive
+ new = [middle; topmiddle; top];
+ len = length(new);
+ oldsteps = linspace(0, 1, len);
+ newsteps = linspace(0, 1, m);
+ newmap = zeros(m, 3);
+
+ for i=1:3
+ % Interpolate over RGB spaces of colormap
+ newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+ end
+
+else
+ % Just negative
+ new = [bottom; botmiddle; middle];
+ len = length(new);
+ oldsteps = linspace(0, 1, len);
+ newsteps = linspace(0, 1, m);
+ newmap = zeros(m, 3);
+
+ for i=1:3
+ % Interpolate over RGB spaces of colormap
+ newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+ end
+
+end
+%
+% m = 64;
+% new = [bottom; botmiddle; middle; topmiddle; top];
+% % x = 1:m;
+%
+% oldsteps = linspace(0, 1, 5);
+% newsteps = linspace(0, 1, m);
+% newmap = zeros(m, 3);
+%
+% for i=1:3
+% % Interpolate over RGB spaces of colormap
+% newmap(:,i) = min(max(interp1(oldsteps, new(:,i), newsteps)', 0), 1);
+% end
+%
+% % set(gcf, 'colormap', newmap), colorbar
\ No newline at end of file
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index a8d4f665..eac0efbf 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -2,7 +2,7 @@ CONTINUE = 1;
JOBNUM = 0; JOBFOUND = 0;
TJOB_SE = []; % Start and end times of jobs
NU_EVOL = []; % evolution of parameter nu between jobs
-MU_EVOL = []; % evolution of parameter nu between jobs
+MU_EVOL = []; % evolution of parameter mu between jobs
ETAB_EVOL= []; %
L_EVOL = []; %
DT_EVOL = []; %
@@ -12,12 +12,16 @@ Ni00_ = []; Ne00_ = [];
GGAMMA_ = [];
PGAMMA_ = [];
PHI_ = [];
+DENS_E_ = [];
+DENS_I_ = [];
+TEMP_E_ = [];
+TEMP_I_ = [];
Ts0D_ = [];
Ts2D_ = [];
Ts5D_ = [];
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);
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
@@ -26,6 +30,10 @@ while(CONTINUE)
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
Pi_new= numel(Pi); Ji_new= numel(Ji);
+ if(Pe_max < Pe_new); Pe_max = Pe_new; end;
+ if(Je_max < Je_new); Je_max = Je_new; end;
+ if(Pi_max < Pi_new); Pi_max = Pi_new; end;
+ if(Ji_max < Ji_new); Ji_max = Ji_new; end;
% If a degree is larger than previous job, put them in a larger array
if (sum([Pe_new, Je_new, Pi_new, Ji_new]>[Pe_old, Je_old, Pi_old, Ji_old]) >= 1)
if W_NAPJ
@@ -44,6 +52,24 @@ while(CONTINUE)
Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
end
+ % If a degree is smaller than previous job, put zero to add. deg.
+ elseif (sum([Pe_new, Je_new, Pi_new, Ji_new]<[Pe_old, Je_old, Pi_old, Ji_old]) >= 1 && Pe_old ~= 1e9)
+ if W_NAPJ
+ tmp = Nipj; sz = size(tmp);
+ Nipj = zeros(cat(1,[Pi_max,Ji_max]',sz(3:end)')');
+ Nipj(1:Pi_new,1:Ji_new,:,:,:) = tmp;
+ tmp = Nepj; sz = size(tmp);
+ Nepj = zeros(cat(1,[Pe_max,Je_max]',sz(3:end)')');
+ Nepj(1:Pe_new,1:Je_new,:,:,:) = tmp;
+ end
+ if W_SAPJ
+ tmp = Sipj; sz = size(tmp);
+ Sipj = zeros(cat(1,[Pi_max,Ji_max]',sz(3:end)')');
+ Sipj(1:Pi_new,1:Ji_new,:,:,:) = tmp;
+ tmp = Sepj; sz = size(tmp);
+ Sepj = zeros(cat(1,[Pe_max,Je_max]',sz(3:end)')');
+ Sepj(1:Pe_new,1:Je_new,:,:,:) = tmp;
+ end
end
if W_GAMMA
@@ -62,7 +88,14 @@ while(CONTINUE)
Ni00_ = cat(3,Ni00_,Ni00);
Ne00_ = cat(3,Ne00_,Ne00);
end
-
+ if W_DENS
+ DENS_E_ = cat(3,DENS_E_,DENS_E);
+ DENS_I_ = cat(3,DENS_I_,DENS_I);
+ end
+ if W_TEMP
+ TEMP_E_ = cat(3,TEMP_E_,TEMP_E);
+ TEMP_I_ = cat(3,TEMP_I_,TEMP_I);
+ end
if W_NAPJ || W_SAPJ
Ts5D_ = cat(1,Ts5D_,Ts5D);
end
@@ -97,6 +130,7 @@ end
GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_;
Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_;
+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_
Sipj = Sipj_; Sepj = Sepj_;
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index e1b57af4..ca3e37b7 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -13,7 +13,7 @@ else
% Setup figure frame
fig = figure('Color','white','Position', [100, 100, 400, 400]);
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;
@@ -27,7 +27,7 @@ fig = figure('Color','white','Position', [100, 100, 400, 400]);
shading interp;
end
set(pclr, 'edgecolor','none'); axis square;
-% caxis([-max(max(abs(FIELD(:,:,n)))),max(max(abs(FIELD(:,:,n))))]);
+ caxis([-1,1]);
xlabel(XNAME); ylabel(YNAME); %colorbar;
title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
,', scaling = ',sprintf('%.1e',scale)]);
diff --git a/matlab/create_mov.m b/matlab/create_mov.m
index 93bc14cb..cebdb625 100644
--- a/matlab/create_mov.m
+++ b/matlab/create_mov.m
@@ -14,7 +14,8 @@ else
% Setup figure frame
figure('Color','white','Position', [100, 100, 400, 400]);
pcolor(X,Y,FIELD(:,:,1)); % to set up
- colormap gray
+% colormap gray
+ colormap(bluewhitered)
axis tight manual % this ensures that getframe() returns a consistent size
if 1
shading interp;
diff --git a/matlab/kernel.m b/matlab/kernel.m
index f0530feb..d0ccac8b 100644
--- a/matlab/kernel.m
+++ b/matlab/kernel.m
@@ -12,4 +12,4 @@ else
kernel_ =0;
end
-end
+end
\ No newline at end of file
diff --git a/matlab/load_params.m b/matlab/load_params.m
index b7fd6315..c74d7b12 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -13,7 +13,8 @@ NZ = h5readatt(filename,'/data/input','nz');
L = h5readatt(filename,'/data/input','Lr');
CLOS = h5readatt(filename,'/data/input','CLOS');
DT_SIM = h5readatt(filename,'/data/input','dt');
-MU = str2num(filename(end-18:end-14));
+% 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';
@@ -26,13 +27,13 @@ else
NON_LIN = 0;
end
-if (CO == -3); CONAME = 'FCDK';
+if (CO == -3); CONAME = 'PADK';
elseif(CO == -2); CONAME = 'SGDK';
elseif(CO == -1); CONAME = 'DGDK';
elseif(CO == 0); CONAME = 'LB';
elseif(CO == 1); CONAME = 'DGGK';
elseif(CO == 2); CONAME = 'SGGK';
-elseif(CO == 3); CONAME = 'FCGK';
+elseif(CO == 3); CONAME = 'PAGK';
end
if (CLOS == 0); CLOSNAME = 'Trunc.';
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
deleted file mode 100644
index 1d3d0db4..00000000
--- a/matlab/photomaton.m
+++ /dev/null
@@ -1,110 +0,0 @@
-if 0
-%% Photomaton : real space
-FIELD = ni00; FNAME = 'ni'; FIELDLTX = '$N_i^{00}$'; XX = RR; YY = ZZ;
-% FIELD = ne00; FNAME = 'ne'; FIELDLTX = '$N_e^{00}$'; XX = RR; YY = ZZ;
-% FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
-% FIELD = dr2phi; FNAME = 'dr2phi'; FIELDLTX = '$\partial_r^2\phi$'; XX = RR; YY = ZZ;
-tf = 10; [~,it1] = min(abs(Ts2D-tf));
-tf = 40; [~,it2] = min(abs(Ts2D-tf));
-tf = 80; [~,it3] = min(abs(Ts2D-tf));
-tf = 100; [~,it4] = min(abs(Ts2D-tf));
-fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
-plt = @(x) x;%./max(max(x));
- subplot(141)
- DATA = plt(FIELD(:,:,it1));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
- subplot(142)
- DATA = plt(FIELD(:,:,it2));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
- subplot(143)
- DATA = plt(FIELD(:,:,it3));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
- subplot(144)
- DATA = plt(FIELD(:,:,it4));
- pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
- suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
- ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
-save_figure
-end
-
-if 0
-%% Photomaton : real space
-% FIELD = Ni00; FNAME = 'fni00'; FIELDLTX = '$\tilde N_i^{00}$'; XX = RR; YY = ZZ;
-% FIELD = Ne00; FNAME = 'fne00'; FIELDLTX = '$\tilde N_e^{00}$'; XX = RR; YY = ZZ;
-FIELD = PHI; FNAME = 'fphi'; FIELDLTX = '$\tilde \phi$'; XX = RR; YY = ZZ;
-tf = 700; [~,it1] = min(abs(Ts2D-tf));
-tf = 800; [~,it2] = min(abs(Ts2D-tf));
-tf = 900; [~,it3] = min(abs(Ts2D-tf));
-tf = 1500; [~,it4] = min(abs(Ts2D-tf));
-XLIM = [0,1.0]; YLIM = [-.5;.5];
-fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 400])
-plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
- subplot(141)
- DATA = plt(FIELD(:,:,it1));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); xlim(XLIM); ylim(YLIM);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
- subplot(142)
- DATA = plt(FIELD(:,:,it2));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); set(gca,'ytick',[]); xlim(XLIM); ylim(YLIM);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
- subplot(143)
- DATA = plt(FIELD(:,:,it3));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$');set(gca,'ytick',[]); xlim(XLIM); ylim(YLIM);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
- subplot(144)
- DATA = plt(FIELD(:,:,it4));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap gray
- xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); set(gca,'ytick',[]); xlim(XLIM); ylim(YLIM);
- title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
- suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
- ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
-save_figure
-end
-
-
-%%
-if 0
-%% Show frame in kspace
-tf = 800; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
-fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600])
-CLIM = [0,1];
- subplot(223); plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none');
- caxis(CLIM); colormap hot
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
- subplot(222);
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none');
- caxis(CLIM); colormap hot
- xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
- subplot(221);
- pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it2))); set(pclr, 'edgecolor','none');
- caxis(CLIM); colormap hot
- xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_e^{00}|$');
- subplot(224);
- colorbar;
- caxis(CLIM); colormap hot
-% pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
-% xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
-save_figure
-end
\ No newline at end of file
diff --git a/matlab/plot_kernels.m b/matlab/plot_kernels.m
new file mode 100644
index 00000000..f7729846
--- /dev/null
+++ b/matlab/plot_kernels.m
@@ -0,0 +1,46 @@
+%% Kernels
+kmax=7;
+nmax=6;
+kr_ = linspace(0,kmax,100);
+
+
+figure
+for n_ = 0:nmax
+ plot(kr_,kernel(n_,kr_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
+end
+ylim_ = ylim;
+plot(kr_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+plot(kr_,J0,'-r','DisplayName','$J_0$');
+legend('show')
+
+%% Bessels and approx
+vperp = linspace(0,1.5,4);
+nmax1=5;
+nmax2=10;
+kmax=7;
+figure
+for i = 1:4
+subplot(2,2,i)
+ v_ = vperp(i);
+ kr_ = linspace(0,kmax,100);
+
+ J0 = besselj(0,kr_*v_);
+ A1 = 1 - kr_.^2*v_^2/4;
+ K1 = zeros(size(kr_));
+ K2 = zeros(size(kr_));
+ for n_ = 0:nmax1
+ K1 = K1 + kernel(n_,kr_).*polyval(LaguerrePoly(n_),v_^2);
+ end
+ for n_ = 0:nmax2
+ K2 = K2 + kernel(n_,kr_).*polyval(LaguerrePoly(n_),v_^2);
+ end
+ plot(kr_,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
+ plot(kr_,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
+ plot(kr_,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kr_,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
+ ylim_ = [-0.5, 1.0];
+ plot(kr_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+ ylim(ylim_); xlabel('$k$')
+ legend('show'); grid on; title(['$v = ',num2str(v_),'$'])
+end
+%%
\ No newline at end of file
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 57b84701..879cd72e 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -1,21 +1,21 @@
%% load profiling
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
+filename_ = ['/misc/HeLaZ_outputs',filename(3:end)];
+CPUTIME = double(h5readatt(filename_,'/data/input','cpu_time'));
+DT_SIM = h5readatt(filename_,'/data/input','dt');
-CPUTIME = double(h5readatt(filename,'/data/input','cpu_time'));
-DT_SIM = h5readatt(filename,'/data/input','dt');
-
-rhs_Tc = h5read(filename,'/profiler/Tc_rhs');
-adv_field_Tc = h5read(filename,'/profiler/Tc_adv_field');
-ghost_Tc = h5read(filename,'/profiler/Tc_ghost');
-clos_Tc = h5read(filename,'/profiler/Tc_clos');
-coll_Tc = h5read(filename,'/profiler/Tc_coll');
-poisson_Tc = h5read(filename,'/profiler/Tc_poisson');
-Sapj_Tc = h5read(filename,'/profiler/Tc_Sapj');
-checkfield_Tc= h5read(filename,'/profiler/Tc_checkfield');
-diag_Tc = h5read(filename,'/profiler/Tc_diag');
-step_Tc = h5read(filename,'/profiler/Tc_step');
-Ts0D = h5read(filename,'/profiler/time');
+rhs_Tc = h5read(filename_,'/profiler/Tc_rhs');
+adv_field_Tc = h5read(filename_,'/profiler/Tc_adv_field');
+ghost_Tc = h5read(filename_,'/profiler/Tc_ghost');
+clos_Tc = h5read(filename_,'/profiler/Tc_clos');
+coll_Tc = h5read(filename_,'/profiler/Tc_coll');
+poisson_Tc = h5read(filename_,'/profiler/Tc_poisson');
+Sapj_Tc = h5read(filename_,'/profiler/Tc_Sapj');
+checkfield_Tc= h5read(filename_,'/profiler/Tc_checkfield');
+diag_Tc = h5read(filename_,'/profiler/Tc_diag');
+step_Tc = h5read(filename_,'/profiler/Tc_step');
+Ts0D = h5read(filename_,'/profiler/time');
missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ...
-coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc;
diff --git a/matlab/setup.m b/matlab/setup.m
index 270c8db0..bf6f7c29 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -38,88 +38,30 @@ MODEL.lambdaD = LAMBDAD;
% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
-if (INIT_PHI && INIT_ZF == 0); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
+if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
INITIAL.INIT_ZF = INIT_ZF;
INITIAL.init_background = (INIT_ZF>0)*ZF_AMP;
INITIAL.init_noiselvl = NOISE0;
INITIAL.iseed = 42;
-INITIAL.selfmat_file = '''null''';
-INITIAL.eimat_file = '''null''';
-INITIAL.iemat_file = '''null''';
INITIAL.mat_file = '''null''';
-if (CO == -3) % Write matrice filename for Full Coulomb
- cmat_pmaxe = 25;
- cmat_jmaxe = 18;
- cmat_pmaxi = 25;
- cmat_jmaxi = 18;
- INITIAL.mat_file = ['''../../../iCa/FC_P_25_J_18_N_200_dk_0.05236_MFLR_0.h5'''];
- INITIAL.selfmat_file = ...
- ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_pamaxx_10.h5'''];
- INITIAL.eimat_file = ...
- ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
- INITIAL.iemat_file = ...
- ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
-elseif (CO == -2) % Write matrice filename for Sugama
- cmat_pmaxe = 10;
- cmat_jmaxe = 5;
- cmat_pmaxi = 10;
- cmat_jmaxi = 5;
-% INITIAL.mat_file = ['''../../../iCa/SG_P_10_J_5_N_200_dk_0.05236_MFLR_0.h5'''];
- INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_10_J_5.h5'''];
- INITIAL.selfmat_file = ...
- ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12.h5'''];
- INITIAL.eimat_file = ...
- ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5'''];
- INITIAL.iemat_file = ...
- ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5'''];
-elseif (CO == 2) % Write matrice filename for Sugama GK
- cmat_pmaxe = 10;
- cmat_jmaxe = 5;
- cmat_pmaxi = 10;
- cmat_jmaxi = 5;
-% INITIAL.mat_file = ['''../../../iCa/SG_P_10_J_5_N_200_dk_0.05236_MFLR_0.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_4_J_2_N_10_kpm_1.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_4_J_2_N_150_kpm_8.h5'''];
+if (abs(CO) == 2) %Sugama operator
INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_10_J_5.h5'''];
- INITIAL.selfmat_file = ...
- ['''../../../iCa/self_Coll_GKE_1_GKI_1_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12_'''];
- INITIAL.eimat_file = ...
- ['''../../../iCa/ei_Coll_GKE_1_GKI_1_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_'''];
- INITIAL.iemat_file = ...
- ['''../../../iCa/ie_Coll_GKE_1_GKI_1_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
- '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
- num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_'''];
-elseif (CO == 3) % Full Coulomb GK
+elseif (abs(CO) == 3) %pitch angle operator
+ INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''];
+elseif (CO == 4) % Full Coulomb GK
disp('Warning, FCGK not implemented yet')
elseif (CO == -1) % DGDK
disp('Warning, DGDK not debugged')
end
% Naming and creating input file
-if (CO == -3); CONAME = 'FCDK';
+if (CO == -3); CONAME = 'PADK';
elseif(CO == -2); CONAME = 'SGDK';
elseif(CO == -1); CONAME = 'DGDK';
elseif(CO == 0); CONAME = 'LB';
elseif(CO == 1); CONAME = 'DGGK';
elseif(CO == 2); CONAME = 'SGGK';
-elseif(CO == 3); CONAME = 'FCGK';
+elseif(CO == 3); CONAME = 'PAGK';
end
if (CLOS == 0); CLOSNAME = 'Trunc.';
elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
@@ -133,7 +75,7 @@ else
end
degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
-
+
degngrad = sprintf(degngrad,[NU,MU]);
if ~NON_LIN; degngrad = ['lin_',degngrad]; end
resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
@@ -141,7 +83,8 @@ gridname = ['L_',num2str(L),'_'];
if (exist('PREFIX','var') == 0); PREFIX = []; end;
if (exist('SUFFIX','var') == 0); SUFFIX = []; end;
PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX];
-BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
+BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
+BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/'];
BASIC.PARAMS = PARAMS;
BASIC.SIMID = SIMID;
BASIC.nrun = 1e8;
@@ -175,6 +118,9 @@ end
if ~exist(BASIC.RESDIR, 'dir')
mkdir(BASIC.RESDIR)
end
+if ~exist(BASIC.MISCDIR, 'dir')
+mkdir(BASIC.MISCDIR)
+end
%% Compile and WRITE input file
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 97dc7b07..19ccdf5d 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -70,9 +70,6 @@ fprintf(fid,[' INIT_ZF =', num2str(INITIAL.INIT_ZF),'\n']);
fprintf(fid,[' init_background =', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed =', num2str(INITIAL.iseed),'\n']);
-fprintf(fid,[' selfmat_file =', INITIAL.selfmat_file,'\n']);
-fprintf(fid,[' eimat_file =', INITIAL.eimat_file,'\n']);
-fprintf(fid,[' iemat_file =', INITIAL.iemat_file,'\n']);
fprintf(fid,[' mat_file =', INITIAL.mat_file,'\n']);
fprintf(fid,'/\n');
diff --git a/wk/ZF_fourier_analysis.m b/wk/ZF_fourier_analysis.m
index 3610ba52..7c37ecad 100644
--- a/wk/ZF_fourier_analysis.m
+++ b/wk/ZF_fourier_analysis.m
@@ -32,30 +32,36 @@ set(gcf, 'Position', [100, 100, 800, 400])
[n,~] = size(Yy);
Pot(:,it) = Yy .* conj(Yy) / n;
end
- [amax, ikmax] = max(mean(Pot,2));
+ [amax, ikZF] = max(mean(Pot,2));
% pclr = pcolor(NN(1:nmax,:),TT(1:nmax,:),Pot(1:nmax,:)); set(pclr, 'edgecolor','none'); hold on;
plot(0:nmax,mean(Pot(1:nmax+1,:),2)/amax,'DisplayName','$\langle\partial_r\phi\rangle_z (k_r)$'); hold on;
- plot([ikmax-1,ikmax-1],[0,1],'--k', 'DisplayName',['$L_z=',num2str(2*pi/kr(ikmax)),'\rho_s$']);
+ plot([ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_z=',num2str(2*pi/kr(ikZF)),'\rho_s$']);
grid on; box on;
title('ZF spatial spectrum')
xlabel('radial mode number'); yticks([]); legend('show')
save_figure
-%% Shear and phi amp phase space
+%% Pred-Pray phase space (A Zonal Flow review, Diamond 2005, Fig 15, Kobayashi 2015)
+
+E_turb = zeros(1,Ns2D); % Time evol. of the turbulence energy (Pred in Kobayashi 2015)
+E_ZF = zeros(1,Ns2D); % Time evol. of the ZF energy (Pray in Kobayashi 2015)
+for it = 1:numel(Ts2D)
+ E_turb(it) = sum(sum((1+KR.^2+KZ.^2).*abs(PHI(:,:,it)).^2))- sum((1+kr.^2).*abs(PHI(:,1,it)).^2);
+ E_ZF(it) = kr(ikZF)^2*abs(PHI(ikZF,1,it)).^2;
+end
fig = figure; FIGNAME = ['phi_shear_phase_space_',PARAMS];
set(gcf, 'Position', [100, 100, 700, 500])
-t1 = Ts2D(end); t0 = 0;
-[~,its2D] = min(abs(Ts2D-t0)); [~,ite2D] = min(abs(Ts2D-t1));
-scatter(phi_maxr_avgz(its2D:ite2D),shear_maxr_avgz(its2D:ite2D),35,Ts2D(its2D:ite2D),'.',...
+scatter(E_ZF*SCALE,E_turb*SCALE,35,Ts2D,'.',...
'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
hold on
-xlabel('$\langle \phi \rangle_z^r$'); ylabel('$\langle dV_E/dr \rangle_z^r$')
+% xlabel('$\langle \phi \rangle_z^r$'); ylabel('$\langle dV_E/dr \rangle_z^r$')
+xlabel('$E_v$'); ylabel('$N$')
grid on; title('ES pot. vs Shear phase space')
% plot(phi_avgr_maxz(its2D:ite2D),shear_avgr_maxz(its2D:ite2D),'-')
% plot(phi_maxr_maxz(its2D:ite2D),shear_maxr_maxz(its2D:ite2D),'-')
% plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
save_figure
-
+clear x_ y_
if 0
%% density and phi phase space
fig = figure; FIGNAME = ['phi_ni_phase_space_',PARAMS];
@@ -74,7 +80,7 @@ grid on; title('ES pot. vs Shear phase space')
end
%% Non zonal quantities
PHI_NZ = PHI;
-PHI_NZ(ikmax-1:ikmax+1,:,:) = 0;
+PHI_NZ(ikZF-1:ikZF+1,:,:) = 0;
phi_nz = zeros(Nr,Nz,Ns2D);
for it = 1:numel(Ts2D)
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 6eebb675..e1896cb8 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,21 +1,28 @@
addpath(genpath('../matlab')) % ... add
+for i_ = 1
% for ETA_ =[0.6:0.1:0.9]
%% Load results
-if 1 % Local results
+if 0% Local results
outfile ='';
outfile ='';
outfile ='';
-outfile ='test_diagnostics/100x50_L_50_P_2_J_1_eta_0.5_nu_5e-03_DGGK_CLOS_0_mu_1e-02';
+outfile ='';
+outfile ='kobayashi/100x50_L_50_P_2_J_1_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00';
+% outfile ='kobayashi/100x50_L_50_P_2_J_1_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00';
+% outfile ='v2.7_P_2_J_1/100x50_L_200_P_2_J_1_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_0e+00';
BASIC.RESDIR = ['../results/',outfile,'/'];
- BASIC.MISCDIR = ['../results/',outfile,'/'];
+ BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
+ CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
+ system(CMD);
end
-if 0 % Marconi results
+if 1% Marconi results
+outfile ='';
outfile ='';
outfile ='';
outfile ='';
outfile ='';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/out.txt';
+% outfile = outcl{i_};
% load_marconi(outfile);
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
@@ -55,7 +62,9 @@ disp('- iFFT')
% IFFT (Lower case = real space, upper case = frequency space)
ne00 = zeros(Nr,Nz,Ns2D); % Gyrocenter density
ni00 = zeros(Nr,Nz,Ns2D);
-dzni00 = zeros(Nr,Nz,Ns2D);
+dzTe = zeros(Nr,Nz,Ns2D);
+dzTi = zeros(Nr,Nz,Ns2D);
+dzni = zeros(Nr,Nz,Ns2D);
np_i = zeros(Nr,Nz,Ns5D); % Ion particle density
si00 = zeros(Nr,Nz,Ns5D);
phi = zeros(Nr,Nz,Ns2D);
@@ -69,19 +78,23 @@ dr2phi = zeros(Nr,Nz,Ns2D);
for it = 1:numel(Ts2D)
NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
- DENS_E_ = DENS_E(:,:,it); DENS_I_ = DENS_I(:,:,it);
- TEMP_E_ = TEMP_E(:,:,it); TEMP_I_ = TEMP_I(:,:,it);
ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz)));
ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz)));
- dzni00(:,:,it) = real(fftshift(ifft2(1i*KZ.*(NI_),Nr,Nz)));
phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz)));
+ drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
+ dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
+ dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
+ if(W_DENS && W_TEMP)
+ DENS_E_ = DENS_E(:,:,it); DENS_I_ = DENS_I(:,:,it);
+ TEMP_E_ = TEMP_E(:,:,it); TEMP_I_ = TEMP_I(:,:,it);
+ dzni(:,:,it) = real(fftshift(ifft2(1i*KZ.*(DENS_I_),Nr,Nz)));
+ dzTe(:,:,it) = real(fftshift(ifft2(1i*KZ.*(TEMP_E_),Nr,Nz)));
+ dzTi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(TEMP_I_),Nr,Nz)));
dens_e (:,:,it) = real(fftshift(ifft2((DENS_E_),Nr,Nz)));
dens_i (:,:,it) = real(fftshift(ifft2((DENS_I_),Nr,Nz)));
temp_e (:,:,it) = real(fftshift(ifft2((TEMP_E_),Nr,Nz)));
temp_i (:,:,it) = real(fftshift(ifft2((TEMP_I_),Nr,Nz)));
- drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
- dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
- dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
+ end
end
% Building a version of phi only 5D sampling times
@@ -118,9 +131,11 @@ PFlux_ri = zeros(1,Ns5D); % Particle flux
% gyrocenter and particle flux from fourier coefficients
GFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Ni00.*conj(PHI),1),2)))*(2*pi/Nr/Nz)^2;
PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz)^2;
+% Heat flux
+Q_RI = -squeeze(mean(mean(dzphi.*temp_i,1),2))';
% Hermite energy spectrum
-epsilon_e_pj = zeros(Npe,Nje,Ns5D);
-epsilon_i_pj = zeros(Npi,Nji,Ns5D);
+epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D);
+epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D);
phi_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
phi_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
@@ -132,11 +147,22 @@ shear_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of shear
shear_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
shear_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
-Ne_norm = zeros(Npe,Nje,Ns5D); % Time evol. of the norm of Napj
-Ni_norm = zeros(Npi,Nji,Ns5D); % .
+Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj
+Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % .
Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2;
-
+% Kperp spectrum interpolation
+%full kperp points
+kperp = reshape(sqrt(KR.^2+KZ.^2),[numel(KR),1]);
+% interpolated kperps
+nk_noAA = floor(2/3*numel(kr));
+kp_ip = kr;
+[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
+[xn,yn] = pol2cart(thg,rg);
+[kz_s, sortIdx] = sort(kz);
+[xc,yc] = meshgrid(kz_s,kr);
+phi_kp_t = zeros(numel(kp_ip),Ns2D);
+%
for it = 1:numel(Ts2D) % Loop over 2D arrays
NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
phi_maxr_maxz(it) = max( max(squeeze(phi(:,:,it))));
@@ -153,8 +179,11 @@ for it = 1:numel(Ts2D) % Loop over 2D arrays
GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
+
+ Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,it))).^2,3)),xn,yn);
+ phi_kp_t(:,it) = mean(Z_rth,2);
end
-
+%
for it = 1:numel(Ts5D) % Loop over 5D arrays
[~, it2D] = min(abs(Ts2D-Ts5D(it)));
Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
@@ -165,7 +194,7 @@ for it = 1:numel(Ts5D) % Loop over 5D arrays
PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lr/Lz;
end
-%% Compute growth rate
+%% Compute primary instability growth rate
disp('- growth rate')
% Find max value of transport (end of linear mode)
[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nr/Nz)^2);
@@ -174,15 +203,33 @@ tstart = 0.1 * Ts2D(itmax); tend = 0.5 * Ts2D(itmax);
[~,its2D_lin] = min(abs(Ts2D-tstart));
[~,ite2D_lin] = min(abs(Ts2D-tend));
-g_ = zeros(Nkr,Nkz);
+g_I = zeros(Nkr,Nkz);
for ikr = 1:Nkr
for ikz = 1:Nkz
- [g_(ikr,ikz), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikr,ikz,its2D_lin:ite2D_lin))));
+ [g_I(ikr,ikz), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikr,ikz,its2D_lin:ite2D_lin))));
+ end
+end
+[gmax_I,ikmax_I] = max(g_I(1,:));
+kmax_I = abs(kz(ikmax_I));
+Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2;
+
+%% Compute secondary instability growth rate
+disp('- growth rate')
+% Find max value of transport (end of linear mode)
+% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nr/Nz)^2);
+% [~,itmax] = min(abs(Ts2D-tmax));
+% tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
+[~,its2D_lin] = min(abs(Ts2D-tstart));
+[~,ite2D_lin] = min(abs(Ts2D-tend));
+
+g_II = zeros(Nkr,Nkz);
+for ikr = 1:Nkr
+ for ikz = 1
+ [g_II(ikr,ikz), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikr,ikz,its2D_lin:ite2D_lin))));
end
end
-[gmax,ikzmax] = max(g_(1,:));
-kzmax = abs(kz(ikzmax));
-Bohm_transport = ETAB/ETAN*gmax/kzmax^2;
+[gmax_II,ikmax_II] = max(g_II(1,:));
+kmax_II = abs(kr(ikmax_II));
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
@@ -198,10 +245,10 @@ subplot(111);
', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
subplot(421);
- for ip = 1:Npe
- for ij = 1:Nje
+ for ip = 1:Pe_max
+ for ij = 1:Je_max
plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$'];
+ plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,...
@@ -210,10 +257,10 @@ subplot(111);
end
grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
subplot(423)
- for ip = 1:Npi
- for ij = 1:Nji
+ for ip = 1:Pi_max
+ for ij = 1:Ji_max
plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$'];
+ plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,...
@@ -229,12 +276,13 @@ subplot(111);
legend(['Gyro. flux';'Part. flux']);
grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
% ylim([0,2.0]);
- if(1)
+ if(0)
subplot(223)
- plot(kz,g_(1,:),'-','DisplayName','Linear growth rate'); hold on;
+ plot(kz,g_I(1,:),'-','DisplayName','Primar. instability'); hold on;
+ plot(kr,g_II(:,1),'x-','DisplayName','Second. instability'); hold on;
plot([max(kz)*2/3,max(kz)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
- grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
- ylim([0,max(g_(1,:))]); xlim([0,max(kz)]);
+ grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
+ ylim([0,max(g_I(1,:))]); xlim([0,max(kz)]);
shearplot = 426; phiplot = 428;
else
shearplot = 223; phiplot = 224;
@@ -260,7 +308,7 @@ end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG = 500; % Averaging time duration
+TAVG = 1000; % Averaging time duration
%Compute steady radial transport
tend = Ts0D(end); tstart = tend - TAVG;
[~,its0D] = min(abs(Ts0D-tstart));
@@ -274,26 +322,30 @@ tend = Ts2D(end); tstart = tend - TAVG;
[~,ite2D] = min(abs(Ts2D-tend));
shear_infty_avg = mean(shear_maxr_avgz(its2D:ite2D));
shear_infty_std = std (shear_maxr_avgz(its2D:ite2D));
+Q_infty_avg = mean(Q_RI(its2D:ite2D))*SCALE;
+Q_infty_std = std(Q_RI(its2D:ite2D))*SCALE;
% plots
fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
- plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
+% yyaxis left
+ plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
plot(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('Transport')
- ylim([0,gamma_infty_avg*5.0]); xlim([0,Ts0D(end)]);
+ grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$')
+ ylim([0,5*abs(gamma_infty_avg)]); xlim([0,Ts0D(end)]);
title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
- legend('show');
-
+% yyaxis right
+% plot(Ts2D,Q_RI*SCALE,'.','DisplayName','$\langle T_i d\phi/dz \rangle_z$'); hold on;
+% ylim([0,5*Q_infty_avg]); xlim([0,Ts0D(end)]); ylabel('$Q_r$')
+% plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Q_infty_avg, '--k',...
+% 'DisplayName',['$Q^{\infty} = $',num2str(Q_infty_avg),'$\pm$',num2str(Q_infty_std)]);
+% legend('show','Location','west')
+ %
subplot(312)
clr = line_colors(1,:);
lstyle = line_styles(1);
-% plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on;
-% plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle \phi\rangle_z$'); hold on;
-% plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle \phi\rangle_r$'); hold on;
-% plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle \phi\rangle_{r,z}$'); hold on;
plot(Ts2D,shear_maxr_maxz,'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
plot(Ts2D,shear_maxr_avgz,'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
plot(Ts2D,shear_avgr_maxz,'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
@@ -301,12 +353,12 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
plot(Ts2D(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([0,Ts0D(end)]);
- grid on; ylabel('Shear amp.'); legend('show');set(gca,'xticklabel',[])
+ grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
subplot(313)
[TY,TX] = meshgrid(r,Ts2D);
% pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_r\phi\rangle_z$') %colorbar;
pclr = pcolor(TX,TY,squeeze(mean(dr2phi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
- caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$r/\rho_s$'); colormap hot
+ caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$r/\rho_s$'); colormap(bluewhitered(256))
save_figure
end
@@ -319,67 +371,94 @@ trange = itstart:itend;
[TY,TX] = meshgrid(r,Ts2D(trange));
fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(211)
+% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
shading interp
colormap hot;
- caxis([0.0,0.8*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
-% caxis([0.0,0.1]);
+ caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
+ caxis([0.0,0.05]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
xticks([]); ylabel('$r/\rho_s$')
- legend('Radial part. transport $\langle n_i^{00}\partial_z\phi\rangle_z$')
+% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
+% subplot(312)
+% pclr = pcolor(TX,TY,squeeze(mean(temp_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+% shading interp
+% % colormap(bluewhitered(256));
+% xticks([]); ylabel('$r/\rho_s$')
+% legend('Radial part. transport $\langle T_i\partial_z\phi\rangle_z$')
+% title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+% ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+% ' $\mu_{hd}=$',num2str(MU)]);
subplot(212)
pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
fieldmax = max(max(mean(abs(drphi(:,:,its2D:ite2D)),2)));
- caxis([-fieldmax,fieldmax]);
+ caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
- legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
+% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
save_figure
end
if 0
%% Averaged shear and Reynold stress profiles
+trange = its2D:ite2D;
+% trange = 100:200;
figure;
-plt = @(x) squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))));
+plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3))));
plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on;
-plot(r,plt(-drphi.*dzphi),'--','DisplayName','$\Pi_\phi$'); hold on;
-plot(r,plt(-drphi.*dzni00),'--','DisplayName','$\Pi_{ni00}$'); hold on;
-plot(r,plt(-drphi.*dzphi-drphi.*dzni00),'--','DisplayName','$\Pi_\phi+\Pi_{ni00}$'); hold on;
-xlim([-60,60]); xlabel('$r/\rho_s$');
+plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on;
+% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on;
+plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on;
+plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on;
+% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on;
+xlim([-L/2,L/2]); xlabel('$r/\rho_s$'); grid on; legend('show')
end
-if 0
+if 1
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 500; tend = 5000;
+tstart = 4000; tend = 4000;
[~,itstart] = min(abs(Ts2D-tstart));
[~,itend] = min(abs(Ts2D-tend));
trange = itstart:itend;
%full kperp points
phi_k_2 = reshape(mean((abs(PHI(:,:,trange))).^2,3),[numel(KR),1]);
kperp = reshape(sqrt(KR.^2+KZ.^2),[numel(KR),1]);
-%interpolated kperps
+% interpolated kperps
nk_noAA = floor(2/3*numel(kr));
-kp_ip = kr(1:nk_noAA);
-[thg, rg] = meshgrid(linspace(-pi/2,pi/2,2*nk_noAA),kp_ip);
+kp_ip = kr;
+[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
[xn,yn] = pol2cart(thg,rg);
-[xc,yc] = meshgrid(sort(kz),kr);
-Z_rth = interp2(xc,yc,mean((abs(PHI(:,:,trange))).^2,3),xn,yn);
-z_k = mean(Z_rth,2);
+[kz_s, sortIdx] = sort(kz);
+[xc,yc] = meshgrid(kz_s,kr);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,trange))).^2,3)),xn,yn);
+phi_kp = mean(Z_rth,2);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(Ni00(:,sortIdx,trange))).^2,3)),xn,yn);
+ni00_kp = mean(Z_rth,2);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(Ne00(:,sortIdx,trange))).^2,3)),xn,yn);
+ne00_kp = mean(Z_rth,2);
%for theorical trends
-kp_1d = linspace(0.1,5,100);
-figure;
-scatter(kperp,phi_k_2,'.'); hold on; grid on;
-plot(kp_ip,z_k);
-plot(kp_1d,1e4*kp_1d.^(-13/3));
-plot(kp_1d,1e4*kp_1d.^(-3)./(1+kp_1d.^2).^2);
-set(gca, 'XScale', 'log');set(gca, 'YScale', 'log');
-xlabel('$k_\perp \rho_s$'); ylabel('$|\phi_k|^2$');
+a1 = phi_kp(2)*kp_ip(2).^(13/3);
+a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
+fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
+% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
+plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
+plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
+plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
+set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on
+xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest')
+title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+', $L=',num2str(L),'$, $N=',num2str(Nr),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+' $\mu_{hd}=$',num2str(MU)]});
xlim([0.1,10]);
+save_figure
+clear Z_rth phi_kp ni_kp Ti_kp
end
%%
-t0 = 000;
+t0 =000;
[~, it02D] = min(abs(Ts2D-t0));
[~, it05D] = min(abs(Ts5D-t0));
skip_ = 1;
@@ -388,20 +467,36 @@ FRAMES_2D = it02D:skip_:numel(Ts2D);
FRAMES_5D = it05D:skip_:numel(Ts5D);
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
+%% part density electron
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+FIELD = real(dens_e); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+create_gif
+% create_mov
+end
+if 0
+%% part temperature electron
+GIFNAME = ['Te',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
+FIELD = real(temp_e); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$T_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
+create_gif
+% create_mov
+end
+if 0
%% part density ion
GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
FIELD = real(dens_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-% create_gif
-create_mov
+create_gif
+% create_mov
end
if 0
%% part temperature ion
GIFNAME = ['Ti',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
FIELD = real(temp_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$T_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-% create_gif
-create_mov
+create_gif
+% create_mov
end
if 0
%% GC Density ion
@@ -425,19 +520,13 @@ GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
+% create_mov
end
if 0
-%% radial particle transport
-GIFNAME = ['gamma_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
-FIELD = real(ni00.*dzphi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$\Gamma_r$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-end
-if 0
-%% Phi fourier
-GIFNAME = ['FFT_phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
-FIELD = ifftshift((abs(PHI)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$|\tilde\phi|$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
+%% shear
+GIFNAME = ['shear_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
+FIELD = -dr2phi; X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
+FIELDNAME = '$s$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
@@ -465,11 +554,11 @@ FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
create_gif_1D_phi
end
if 0
-%% phi @ kz = 0
-GIFNAME = ['phi_kz0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
-FIELD =squeeze(log10(abs(PHI(:,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
-X = kr; T = Ts2D; YMIN = -0; YMAX = 6; XMIN = min(kr); XMAX = max(kr);
-FIELDNAME = '$\log_{10}|\tilde\phi(k_z=0)|$'; XNAME = '$k_r\rho_s$';
+%% phi kperp spectrum
+GIFNAME = ['phi2_kp',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
+FIELD =log10(phi_kp_t); linestyle = '-'; FRAMES = FRAMES_2D;
+X = kp_ip; T = Ts2D; YMIN = -20; YMAX = 10; XMIN = min(kr); XMAX = max(kr);
+FIELDNAME = '$\log_{10}|\tilde\phi_k|^2$'; XNAME = '$k_r\rho_s$';
create_gif_1D
end
if 0
@@ -512,4 +601,4 @@ create_gif_5D
end
%%
ZF_fourier_analysis
-% end
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index a52c803f..c2b825b0 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,12 +1,9 @@
%% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_60_P_6_J_3_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_1e-03/out.txt')
%% Functions to modify preexisting fort.90 input file and launch on marconi
function [] = continue_run(outfilename)
- EXECNAME = 'helaz_2.6';
+ EXECNAME = 'helaz_2.73';
%% CLUSTER PARAMETERS
CLUSTER.PART = 'prod'; % dbg or prod
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
@@ -44,15 +41,25 @@ function [] = continue_run(outfilename)
A{5} = sprintf(' RESTART = .true.');
J2L = 0;
else
- line = A{33};
+ line = A{35};
line = line(end-2:end);
if(line(1) == '='); line = line(end); end;
- J2L = str2num(line) + 1;
+ J2L = str2num(line);
end
% Change job 2 load in fort.90
- A{33} = [' job2load = ',num2str(J2L)];
- disp(A{33})
-
+ A{35} = [' job2load = ',num2str(0)];
+ disp(A{35})
+ % Change time step
+ A{3} = [' dt = 0.01'];
+ % Increase endtime
+ A{4} = [' tmax = 10000'];
+ % Put non linear term back
+ A{41} = [' NL_CLOS = -1'];
+ % change HD
+ line_= A{43};
+ mu_old = str2num(line_(13:end));
+ A{43} = [' mu = ',num2str(1e-3)];
+
% Rewrite fort.90
fid = fopen('fort.90', 'w');
for i = 1:numel(A)
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 58855861..5838766d 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,6 +1,6 @@
-for NU = [1e-1]
-for ETAB = [0.5]
-for CO = [1]
+for NU = [0.14]
+for ETAB = [1/1.4]
+for CO = [3]
%clear all;
addpath(genpath('../matlab')) % ... add
default_plots_options
@@ -18,7 +18,7 @@ NU_HYP = 0.0; % Hyperdiffusivity coefficient
LAMBDAD = 0.0;
NOISE0 = 1.0e-5;
%% GRID PARAMETERS
-N = 100; % Frequency gridpoints (Nkr = N/2)
+N = 50; % Frequency gridpoints (Nkr = N/2)
L = 100; % Size of the squared frequency domain
KREQ0 = 1; % put kr = 0
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
@@ -26,20 +26,21 @@ MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARMETERS
TMAX = 100; % Maximal time unit
DT = 1e-2; % Time step
-SPS0D = 0.5; % Sampling per time unit for 2D arrays
+SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/2; % Sampling per time unit for 5D arrays
+SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 00;
%% OPTIONS
-SIMID = 'test'; % Name of the simulation
-% SIMID = 'kobayashi_lin'; % Name of the simulation
+% SIMID = 'v2.7_lin_analysis'; % Name of the simulation
+SIMID = 'kobayashi_2015_fig2'; % Name of the simulation
% SIMID = 'v2.6_lin_analysis'; % Name of the simulation
NON_LIN = 0 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
% CO = 2;
+INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
NL_CLOS = 0; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
@@ -47,11 +48,13 @@ INIT_PHI= 0; % Start simulation with a noisy phi
INIT_ZF = 0; % Start simulation with a noisy phi
%% OUTPUTS
W_DOUBLE = 0;
-W_GAMMA = 1;
-W_PHI = 1;
+W_GAMMA = 0;
+W_PHI = 0;
W_NA00 = 1;
W_NAPJ = 1;
W_SAPJ = 0;
+W_DENS = 0;
+W_TEMP = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -65,10 +68,10 @@ MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
%% PARAMETER SCANS
if 1
%% Parameter scan over PJ
-% PA = [2, 4, 6, 8, 10];
-% JA = [1, 2, 3, 4, 5];
-PA = [4];
-JA = [2];
+PA = [2, 6, 10];
+JA = [1, 3, 5];
+% PA = [2 4];
+% JA = [1 2];
DTA= DT./sqrt(JA)/4;
% DTA= DT;
mup_ = MU_P;
@@ -78,7 +81,7 @@ param_name = 'PJ';
gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
gamma_Nipj = zeros(Nparam,floor(N/2)+1);
Bohm_transport = zeros(Nparam,1);
-Ni00_ST = zeros(Nparam,floor(N/2)+1,SPS2D*TMAX);
+Ni00_ST = zeros(Nparam,floor(N/2)+1,floor(SPS2D*TMAX));
for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
@@ -125,28 +128,15 @@ plt = @(x) x;
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*kr),plt(gamma_Ni00(i,:)),...
+ semilogx(plt(SCALE*kr(2:numel(kr))),plt(gamma_Ni00(i,2:end)),...
'Color',clr,...
- 'LineStyle',linestyle{1},...
- 'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
+ 'LineStyle',linestyle{1},'Marker','^',...
+ 'DisplayName',['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
hold on;
end
grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kr)]);
- title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
- legend('show')
-% subplot(212)
-% 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*kr),plt(gamma_Nipj(i,:)),...
-% 'Color',clr,...
-% 'LineStyle',linestyle{1},...
-% 'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
-% hold on;
-% end
-% grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(\max_{pj}N_i^{pj})\rho_s/c_s$'); xlim([0.0,max(kr)]);
-% title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
-% legend('show')
+ title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
+ legend('show'); xlim([0.01,10])
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
end
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 0c7c81ea..ee1d4021 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,8 +1,5 @@
addpath(genpath('../matlab')) % ... add
%% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.5_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/100x50_L_50_P_10_J_5_eta_0.71429_nu_5e-03_DGGK_CLOS_0_mu_1e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_60_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_1e-03/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_1e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 7d099685..a7a196a5 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,10 +4,10 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 5.0e-3; % Collision frequency
-ETAB = 0.5;%1/1.4; % Magnetic gradient
+NU = 0.0141; % Collision frequency
+ETAB = 1/1.4; % Magnetic gradient
ETAN = 1.0; % Density gradient
-NU_HYP = 1.0;
+NU_HYP = 0.0;
%% GRID PARAMETERS
N = 100; % Frequency gridpoints (Nkr = N/2)
L = 50; % Size of the squared frequency domain
@@ -16,22 +16,22 @@ J = 1;
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 500; % Maximal time unit
-DT = 2e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
+TMAX = 1000; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/50; % Sampling per time unit for 5D arrays
+SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 0;
%% OPTIONS AND NAMING
% Collision operator
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO =1;
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
+CO = 3;
CLOS = 0; % Closure model (0: =0 truncation)
NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'test_diagnostics'; % Name of the simulation
-% SIMID = 'kobayashi'; % Name of the simulation
+% SIMID = 'test_restart'; % Name of the simulation
+SIMID = 'kobayashi'; % Name of the simulation
% SIMID = ['v2.7_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
NON_LIN = 1; % activate non-linearity (is cancelled if KREQ0 = 1)
INIT_ZF = 0; ZF_AMP = 0.0;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index eab3af94..1e8aabdb 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -2,13 +2,14 @@ clear all;
addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_2.71';
+ EXECNAME = 'helaz_2.73';
for ETAB = [0.6]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLUSTER PARAMETERS
-CLUSTER.PART = 'prod'; % dbg or prod
+% CLUSTER.PART = 'prod'; % dbg or prod
+CLUSTER.PART = 'dbg';
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
CLUSTER.MEM = '128GB'; % Memory
@@ -19,35 +20,35 @@ NP_KR = 24; % MPI processes along kr
%% PHYSICAL PARAMETERS
NU = 1e-3; % Collision frequency
% ETAB = 0.7; % Magnetic gradient
-NU_HYP = 1.0; % Hyperdiffusivity coefficient
-NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 1;
+NU_HYP = 0.1; % Hyperdiffusivity coefficient
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
+CO = 2;
INIT_ZF = 0; ZF_AMP = 0.0;
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 120; % Size of the squared frequency domain
-P = 10; % Electron and Ion highest Hermite polynomial degree
-J = 5; % Electron and Ion highest Laguerre polynomial degree
+L = 60; % Size of the squared frequency domain
+P = 6; % Electron and Ion highest Hermite polynomial degree
+J = 3; % Electron and Ion highest Laguerre polynomial degree
MU_P = 0.0;% Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0;% Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 5000; % Maximal time unit
-DT = 1e-5; % Time step
+TMAX = 10000; % Maximal time unit
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1/4; % Sampling per time unit for 2D arrays
-SPS5D = 1/100; % Sampling per time unit for 5D arrays
+SPS5D = 1/300; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 0;
%% Naming
-% SIMID = 'kobayashi'; % Name of the simulation
+SIMID = 'kobayashi'; % Name of the simulation
% SIMID = 'test'; % Name of the simulation
-SIMID = ['v2.7_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
+% SIMID = ['v2.7_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
PREFIX =[];
% PREFIX = sprintf('%d_%d_',NP_P, NP_KR);
%% Options
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_PHI= 1; % Start simulation with a noisy phi and moments
%% OUTPUTS
@@ -57,10 +58,12 @@ W_PHI = 1;
W_NA00 = 1;
W_NAPJ = 1;
W_SAPJ = 0;
+W_DENS = 1;
+W_TEMP = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% fixed parameters (for current study)
-KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
+KR0KH = 0; A0KH = 0; % Background phi mode
KREQ0 = 0; % put kr = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
index 947dcf31..b7ef07f7 100644
--- a/wk/open_figure_script.m
+++ b/wk/open_figure_script.m
@@ -1,28 +1,26 @@
simname_ = '';
fname_='';
%% Marconi output file
-% fname_='';
-% fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
-fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+fname_='';
+fname_='';
+fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
simname_ = fname_(54:end-8);
%%
% simname_ = '';
% simname_ = '';
% simname_ = '';
-% simname_ = 'v2.7_P_4_J_2/100x50_L_50_P_4_J_2_eta_0.7_nu_5e-02_DGGK_CLOS_0_mu_1e-02';
+% simname_ = '';
+% simname_ = '';
+% simname_ = 'v2.7_P_2_J_1/100x50_L_200_P_2_J_1_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_0e+00';
+simname_ = 'v2.7_P_2_J_1/100x50_L_100_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_0e+00';
-% simname_ = 'v2.6_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-01';
-% simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
-% simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
-% simname_ = 'v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
-% simname_ = 'v2.6_P_4_J_2/200x100_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
%%
-% figname_ = '/fig/ZF_transport_drphi_';
+figname_ = '/fig/ZF_transport_drphi_';
% figname_ = '/fig/space_time_';
-figname_ = '/fig/phi_shear_phase_space_';
+% figname_ = '/fig/phi_shear_phase_space_';
path_ = '../results/';
diff --git a/wk/photomaton.m b/wk/photomaton.m
new file mode 100644
index 00000000..32b5454d
--- /dev/null
+++ b/wk/photomaton.m
@@ -0,0 +1,90 @@
+if 0
+%% Photomaton : real space
+% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = '$n_i^{00}$'; XX = RR; YY = ZZ;
+FIELD = ne00; FNAME = 'ne00'; FIELDLTX = '$n_e^{00}$'; XX = RR; YY = ZZ;
+% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = '$n_i$'; XX = RR; YY = ZZ;
+% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = '$n_e$'; XX = RR; YY = ZZ;
+% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = '$T_i$'; XX = RR; YY = ZZ;
+% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = '$T_e$'; XX = RR; YY = ZZ;
+% FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
+% FIELD = drphi; FNAME = 'ZF'; FIELDLTX = '$u^{ZF}_z$'; XX = RR; YY = ZZ;
+% FIELD = -dr2phi; FNAME = 'shear'; FIELDLTX = '$s$'; XX = RR; YY = ZZ;
+TNAME = [];
+tf = 250; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 450; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1200; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1500; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
+plt = @(x) x./max(max(x));
+ subplot(141)
+ DATA = plt(FIELD(:,:,it1));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap(bluewhitered); caxis([-1,1]);
+ xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
+ subplot(142)
+ DATA = plt(FIELD(:,:,it2));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap(bluewhitered); caxis([-1,1]);
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
+ subplot(143)
+ DATA = plt(FIELD(:,:,it3));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap(bluewhitered); caxis([-1,1]);
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
+ subplot(144)
+ DATA = plt(FIELD(:,:,it4));
+ pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ colormap(bluewhitered); caxis([-1,1]);
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
+ suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+ ' $\mu_{hd}=$',num2str(MU)]);
+save_figure
+end
+
+if 0
+%% Photomaton : quiver ExB velocity
+figure
+skip = 2;
+plt = @(x) x./max(max(x));
+FNAME = 'ZF'; FIELDLTX = '$\bm{u}^{ZF}$'; XX = RR; YY = 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(XX,YY,plt(ni00(:,:,it1))); set(pclr, 'edgecolor','none');
+hold on
+quiver((XX(1:skip:end,1:skip:end)),(YY(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
+
+
+%%
+if 0
+%% Show frame in kspace
+tf = 800; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
+fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600])
+CLIM = [0,1];
+ subplot(223); plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
+ subplot(222);
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
+ subplot(221);
+ pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it2))); set(pclr, 'edgecolor','none');
+ caxis(CLIM); colormap hot
+ xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_e^{00}|$');
+ subplot(224);
+ colorbar;
+ caxis(CLIM); colormap hot
+% pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
+% xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
+save_figure
+end
\ No newline at end of file
--
GitLab