From 6ada1c908c4f5ac066a26cf5f02636fc3b35c8d3 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 30 Aug 2021 15:52:38 +0200
Subject: [PATCH] scripts update
---
matlab/compile_results.m | 6 +-
matlab/dnjs.m | 6 +-
matlab/dnjs_fact.m | 16 ++
matlab/dnjs_stir.m | 18 ++
matlab/plot_kernels.m | 204 ++++++++++++++++--
matlab/plots/plot_kperp_spectrum.m | 33 ++-
.../plots/plot_radial_transport_and_shear.m | 10 +-
matlab/post_processing.m | 113 +++++-----
matlab/setup.m | 3 +-
wk/HD_study.m | 120 +++++++++--
wk/analysis_3D.m | 126 +++++++----
wk/continue_multiple_runs_marconi.m | 37 ++--
wk/linear_study.m | 10 +-
wk/load_multiple_outputs_marconi.m | 5 +-
wk/local_run.m | 24 +--
wk/marconi_run.m | 28 +--
wk/plot_cosol_mat.m | 71 +++---
17 files changed, 593 insertions(+), 237 deletions(-)
create mode 100644 matlab/dnjs_fact.m
create mode 100644 matlab/dnjs_stir.m
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index b79637a5..44177cf2 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -3,7 +3,7 @@ JOBNUM = JOBNUMMIN; JOBFOUND = 0;
TJOB_SE = []; % Start and end times of jobs
NU_EVOL = []; % evolution of parameter nu between jobs
MU_EVOL = []; % evolution of parameter mu between jobs
-ETAB_EVOL= []; %
+ETAN_EVOL= []; %
L_EVOL = []; %
DT_EVOL = []; %
% FIELDS
@@ -113,7 +113,7 @@ while(CONTINUE)
TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
NU_EVOL = [NU_EVOL NU NU];
MU_EVOL = [MU_EVOL MU MU];
- ETAB_EVOL = [ETAB_EVOL ETAB ETAB];
+ ETAN_EVOL = [ETAN_EVOL ETAN ETAN];
L_EVOL = [L_EVOL L L];
DT_EVOL = [DT_EVOL DT_SIM DT_SIM];
@@ -135,5 +135,5 @@ clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_
Sipj = Sipj_; Sepj = Sepj_;
clear Sipj_ Sepj_
-JOBNUM = LASTJOB
+JOBNUM = LASTJOB;
filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
\ No newline at end of file
diff --git a/matlab/dnjs.m b/matlab/dnjs.m
index d333c1b2..1253d901 100644
--- a/matlab/dnjs.m
+++ b/matlab/dnjs.m
@@ -1,11 +1,11 @@
function [ result ] = dnjs( n, j, s )
+% Compute the dnjs from Ln*Lj = sum_s dnjs Ls with Laguerre coeffs
% sort in order to compute only once the laguerre coeff
Coeffs = sort([n,j,s]);
% last element of Coeffs is the larger one
L3 = flip(LaguerrePoly(Coeffs(end)));
-% retrive smaller order coeff by taking the firsts (flipped array)
- L2 = L3(1:Coeffs(2)+1);
- L1 = L3(1:Coeffs(1)+1);
+ L2 = flip(LaguerrePoly(Coeffs(end-1)));
+ L1 = flip(LaguerrePoly(Coeffs(end-2)));
% build a factorial array to compute once every needed factorials
Factar = zeros(n+j+s+1,1);
diff --git a/matlab/dnjs_fact.m b/matlab/dnjs_fact.m
new file mode 100644
index 00000000..ee5faf82
--- /dev/null
+++ b/matlab/dnjs_fact.m
@@ -0,0 +1,16 @@
+function [res] = dnjs_fact(n,j,s)
+% Compute the dnjs from Ln*Lj = sum_s dnjs Ls with factorials only
+ binomial = @(n,k) factorial(n)./factorial(n-k)./factorial(k);
+ res = 0;
+ for n1 = 0:n
+ for j1 = 0:j
+ for s1 = 0:s
+ res = res ...
+ +(-1).^(n1+j1+s1)...
+ .* binomial(n,n1) .* binomial(j,j1) .* binomial(s,s1)...
+ ./ factorial(n1) ./ factorial(j1) ./ factorial(s1)...
+ .* factorial(n1+j1+s1);
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/dnjs_stir.m b/matlab/dnjs_stir.m
new file mode 100644
index 00000000..013ac5f3
--- /dev/null
+++ b/matlab/dnjs_stir.m
@@ -0,0 +1,18 @@
+function [res] = dnjs_stir(n,j,s)
+% Compute the dnjs from Ln*Lj = sum_s dnjs Ls with stirling formula
+ stirling = @(n,k) sqrt(2*pi).*n.^(n+.5).*exp(-n) .* (n~=0) + (n==0);
+ bin_stir = @(n,k) stirling(n)./stirling(n-k)./stirling(k);
+
+ res = 0;
+ for n1 = 0:n
+ for j1 = 0:j
+ for s1 = 0:s
+ res = res ...
+ +(-1).^(n1+j1+s1)...
+ .* bin_stir(n,n1) .* bin_stir(j,j1) .* bin_stir(s,s1)...
+ ./ stirling(n1) ./ stirling(j1) ./ stirling(s1)...
+ .* stirling(n1+j1+s1);
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/plot_kernels.m b/matlab/plot_kernels.m
index f7eede8a..400e63cc 100644
--- a/matlab/plot_kernels.m
+++ b/matlab/plot_kernels.m
@@ -1,18 +1,22 @@
+if 0
%% Kernels
kmax=7;
nmax=6;
-kx_ = linspace(0,kmax,100);
+kx = linspace(0,kmax,100);
figure
for n_ = 0:nmax
- plot(kx_,kernel(n_,kx_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
+ plot(kx,kernel(n_,kx),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
end
ylim_ = ylim;
-plot(kx_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
-plot(kx_,J0,'-r','DisplayName','$J_0$');
+plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+plot(kx,J0,'-r','DisplayName','$J_0$');
legend('show')
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if 0
%% Bessels and approx
vperp = linspace(0,1.5,4);
nmax1=5;
@@ -22,25 +26,191 @@ figure
for i = 1:4
subplot(2,2,i)
v_ = vperp(i);
- kx_ = linspace(0,kmax,100);
+ kx = linspace(0,kmax,100);
- J0 = besselj(0,kx_*v_);
- A1 = 1 - kx_.^2*v_^2/4;
- K1 = zeros(size(kx_));
- K2 = zeros(size(kx_));
+ J0 = besselj(0,kx*v_);
+ A1 = 1 - kx.^2*v_^2/4;
+ K1 = zeros(size(kx));
+ K2 = zeros(size(kx));
for n_ = 0:nmax1
- K1 = K1 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2);
+ K1 = K1 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
for n_ = 0:nmax2
- K2 = K2 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2);
+ K2 = K2 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
- plot(kx_,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
- plot(kx_,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
- plot(kx_,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
- plot(kx_,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kx,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
+ plot(kx,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
+ plot(kx,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kx,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
ylim_ = [-0.5, 1.0];
- plot(kx_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+ plot(kx(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
+%%
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if 0
+%% from Sum_n Kn x Ln x Lj to Sum_n Kn Sum_s dnjs x Ls
+vperp = [5];
+kx = linspace(0,10,200);
+Jmax = 15;
+j = 10;
+fig = figure; set(gcf, 'Position', [100, 100, numel(vperp)*300, 250])
+% suptitle(['$J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
+Kn = @(x__,y__) kernel(x__,y__);
+for i = 1:numel(vperp)
+subplot(1,numel(vperp),i)
+ v_ = vperp(i);
+ % J0 x Lj
+ J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
+ % Taylor exp of J0
+ A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
+ % Sum_n^Nmax Kn x Ln x Lj
+ KLL = zeros(size(kx));
+ for n_ = 0:Jmax
+ KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
+ end
+ KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL1 = zeros(size(kx));
+ for n_ = 0:Jmax-j
+ tmp_ = 0;
+ for s_ = 0:n_+j
+ tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL1 = KdL1 + Kn(n_,kx).*tmp_;
+ end
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL2 = zeros(size(kx));
+ for n_ = 0:Jmax
+ tmp_ = 0;
+ for s_ = 0:min(Jmax,n_+j)
+ tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL2 = KdL2 + Kn(n_,kx).*tmp_;
+ end
+
+ % plots
+ plot(kx,J0Lj,'-k','DisplayName','$J_0 L_j$'); hold on;
+ plot(kx,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L^n L_j$']);
+ plot(kx,KdL1,'-.','DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L^s$']);
+ ylim_ = ylim;
+ plot(kx,A1,'-','DisplayName','$(1 - k^2 v^2/4) L_j$'); hold on;
+ plot(kx,KdL2,'--','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L^s$']);
+% plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+ ylim(ylim_);
+ xlabel('$k_{\perp}$')
+ legend('show');
+ grid on; title(['$v = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
+end
+FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(j),'.eps'];
+% saveas(fig,FIGNAME,'epsc');
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if 1
+%% Convergence analysis
+kx = linspace(0,10,200);
+v_A = linspace(0,1,10);
+J_A = 1:15;
+[YY,XX] = meshgrid(v_A,J_A);
+klimLL = zeros(size(XX));
+klimdL1= zeros(size(XX));
+Kn = @(x__,y__) kernel(x__,y__);
+fig = figure; set(gcf, 'Position', [100, 100, 500, 400]);
+for j = 5
+for ij_ = 1:numel(J_A)
+ iv_ = 1;
+ for v_ = v_A
+ eps = abs(0.01*polyval(LaguerrePoly(j),v_^2));
+ Jmax = J_A(ij_);
+ % J0 x Lj
+ J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
+ % Taylor exp of J0
+ A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
+ % Sum_n^Nmax Kn x Ln x Lj
+ KLL = zeros(size(kx));
+ for n_ = 0:Jmax
+ KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
+ end
+ KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL1 = zeros(size(kx));
+ for n_ = 0:Jmax-j
+ tmp_ = 0;
+ for s_ = 0:n_+j
+ tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL1 = KdL1 + Kn(n_,kx).*tmp_;
+ end
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL2 = zeros(size(kx));
+ for n_ = 0:Jmax
+ tmp_ = 0;
+ for s_ = 0:min(Jmax,n_+j)
+ tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL2 = KdL2 + Kn(n_,kx).*tmp_;
+ end
+ % errors
+ err_ = abs(J0Lj-KLL);
+ idx_ = 1;
+ while(err_(idx_)< eps && idx_ < numel(err_))
+ idx_ = idx_ + 1;
+ end
+ klimLL(ij_,iv_) = kx(idx_);
+ err_ = abs(J0Lj-KdL1);
+ idx_ = 1;
+ while(err_(idx_)< eps && idx_ < numel(err_))
+ idx_ = idx_ + 1;
+ end
+ klimdL1(ij_,iv_) = kx(idx_);
+ iv_ = iv_ + 1;
+ end
+end
+% plot
+% plot(JmaxA,klimLL,'o-'); hold on
+% plot(J_A,klimdL1,'o-'); hold on
+end
+surf(XX,YY,klimdL1)
+xlabel('$J_{max}$'); ylabel('$v_{\perp}$');
+title(['$k$ s.t. $\epsilon=1\%$, $j=',num2str(j),'$'])
+ylim([0,1]); grid on;
+end
+if 0
+%% Test Lj Ln = sum_s dnjs Ls
+j = 10;
+n = 10;
+smax = n+j;
+vperp = linspace(0,5,20);
+LjLn = zeros(size(vperp));
+dnjsLs = zeros(size(vperp));
+dnjsLs_bin = zeros(size(vperp));
+dnjsLs_stir = zeros(size(vperp));
+
+i=1;
+for x_ = vperp
+ LjLn(i) = polyval(LaguerrePoly(j),x_)*polyval(LaguerrePoly(n),x_);
+ dnjsLs(i) = 0;
+ dnjsLs_bin(i) = 0;
+ dnjsLs_stir(i) = 0;
+ for s_ = 0:smax
+ dnjsLs(i) = dnjsLs(i) + dnjs(n,j,s_)*polyval(LaguerrePoly(s_),x_);
+ dnjsLs_bin(i) = dnjsLs_bin(i) + dnjs_fact(n,j,s_)*polyval(LaguerrePoly(s_),x_);
+ dnjsLs_stir(i) = dnjsLs_stir(i) + dnjs_stir(n,j,s_)*polyval(LaguerrePoly(s_),x_);
+ end
+ i = i+1;
+end
+
+figure
+plot(vperp,LjLn); hold on;
+plot(vperp,dnjsLs,'d')
+plot(vperp,dnjsLs_bin,'o')
+plot(vperp,dnjsLs_stir/dnjsLs_stir(1),'x')
+% We see that starting arround j = 18, n = 0, the stirling formula is the only
+% approximation that does not diverge from the target function. However it
+% performs badly for non 0 n...
+end
+
+
+
diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plots/plot_kperp_spectrum.m
index 2370dd02..238faa03 100644
--- a/matlab/plots/plot_kperp_spectrum.m
+++ b/matlab/plots/plot_kperp_spectrum.m
@@ -2,7 +2,7 @@
[~,itend] = min(abs(Ts3D-tend));
trange = itstart:itend;
%full kperp points
-phi_k_2 = reshape(mean(mean(abs(PHI(:,:,:,trange)),3),4).^2,[numel(KX),1]);
+phi_k_2 = reshape(mean(mean(abs(FIELD(:,:,:,trange)),3),4).^2,[numel(KX),1]);
kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
% interpolated kperps
nk_noAA = floor(2/3*numel(kx));
@@ -11,27 +11,26 @@ kp_ip = kx;
[xn,yn] = pol2cart(thg,rg);
[ky_s, sortIdx] = sort(ky);
[xc,yc] = meshgrid(ky_s,kx);
-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);
+Z_rth = interp2(xc,yc,squeeze(mean((abs(FIELD(:,sortIdx,trange))).^2,3)),xn,yn);
+field_kp = mean(Z_rth,2);
%for theorical trends
-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])
+a1 = field_kp(2)*kp_ip(2).^(13/3);
+a2 = field_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
+fig = figure; FIGNAME = ['cascade','_',FNAME,'_',PARAMS];set(gcf, 'Position', [100, 100, 800, 300])
% 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,field_kp,'^','DisplayName',['$\langle|',FIELDLTX,'|^2\rangle_{k_\perp}$']); hold on;
+if TRENDS
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')
+end
+if LOGSCALE
+ set(gca, 'XScale', 'log');set(gca, 'YScale', 'log');
+ xlim([0.1,10]);
+end
+grid on
+xlabel('$k_\perp \rho_s$'); legend('show','Location','northeast')
title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-' $\mu_{hd}=$',num2str(MU)]});
-xlim([0.1,10]);
+' $\mu_{hd}=$',num2str(MU),', $',num2str(round(tstart)),'<t<',num2str(round(tend)),'$']});
save_figure
clear Z_rth phi_kp ni_kp Ti_kp
\ No newline at end of file
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index ba865250..734551cc 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -29,16 +29,16 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
lstyle = line_styles(1);
% plt = @(x_) mean(x_,1);
plt = @(x_) x_(1,:);
- plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
- plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
- plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
- plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on;
+ plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(s_\phi)$'); hold on;
+ plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle s_\phi\rangle_y$'); hold on;
+ plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle s_\phi\rangle_x$'); hold on;
+ plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{x,y}$'); hold on;
plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
subplot(313)
[TY,TX] = meshgrid(x,Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
+ pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_x^2\phi\rangle_y$)') %colorbar;
caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
save_figure
\ No newline at end of file
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index ea25ef89..5b50bd39 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -1,4 +1,3 @@
-
%% Retrieving max polynomial degree and sampling info
Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
@@ -15,10 +14,13 @@ Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
KPERP2 = KY.^2+KX.^2;
[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
+[KY_XY,KX_XY] = meshgrid(ky,kx);
+[KZ_XZ,KX_XZ] = meshgrid(z,kx);
+[KZ_YZ,KY_YZ] = meshgrid(z,ky);
Lk = max(Lkx,Lky);
Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z);
-dx = 2*pi/Lk; dy = 2*pi/Lk; dz = q0*2*pi/Nz;
+dx = 2*pi/Lk; dy = 2*pi/Lk; dz = 2*pi/Nz;
x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
z = dz * (1:Nz);
@@ -30,72 +32,64 @@ z = dz * (1:Nz);
disp('Analysis :')
disp('- iFFT')
% IFFT (Lower case = real space, upper case = frequency space)
-ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density
-ni00 = zeros(Nx,Ny,Nz,Ns3D);
-dzTe = zeros(Nx,Ny,Nz,Ns3D);
-dzTi = zeros(Nx,Ny,Nz,Ns3D);
-dzni = zeros(Nx,Ny,Nz,Ns3D);
-np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density
-si00 = zeros(Nx,Ny,Nz,Ns5D);
-phi = zeros(Nx,Ny,Nz,Ns3D);
-dens_e = zeros(Nx,Ny,Nz,Ns3D);
-dens_i = zeros(Nx,Ny,Nz,Ns3D);
-temp_e = zeros(Nx,Ny,Nz,Ns3D);
-temp_i = zeros(Nx,Ny,Nz,Ns3D);
-drphi = zeros(Nx,Ny,Nz,Ns3D);
-dzphi = zeros(Nx,Ny,Nz,Ns3D);
-dr2phi = zeros(Nx,Ny,Nz,Ns3D);
+ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density
+ni00 = zeros(Nx,Ny,Nz,Ns3D); % "
+phi = zeros(Nx,Ny,Nz,Ns3D); % Electrostatic potential
+Z_phi = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
+dens_e = zeros(Nx,Ny,Nz,Ns3D); % Particle density
+dens_i = zeros(Nx,Ny,Nz,Ns3D); %"
+Z_n_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
+Z_n_i = zeros(Nx,Ny,Nz,Ns3D); %"
+temp_e = zeros(Nx,Ny,Nz,Ns3D); % Temperature
+temp_i = zeros(Nx,Ny,Nz,Ns3D); % "
+Z_T_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
+Z_T_i = zeros(Nx,Ny,Nz,Ns3D); %"
+dyTe = zeros(Nx,Ny,Nz,Ns3D); % Various derivatives
+dyTi = zeros(Nx,Ny,Nz,Ns3D); % "
+dyni = zeros(Nx,Ny,Nz,Ns3D); % "
+dxphi = zeros(Nx,Ny,Nz,Ns3D); % "
+dyphi = zeros(Nx,Ny,Nz,Ns3D); % "
+dx2phi = zeros(Nx,Ny,Nz,Ns3D); % "
for it = 1:numel(Ts3D)
for iz = 1:numel(z)
NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
- ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
- ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
- phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
- drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny)));
- dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
- dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
- if(W_DENS && W_TEMP)
+ ne00 (:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
+ ni00 (:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
+ phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
+ Z_phi (:,:,iz,it) = real(fftshift(ifft2((PH_.*(KY==0)),Nx,Ny)));
+ dxphi (:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny)));
+ dx2phi(:,:,iz,it) = real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
+ dyphi (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
+ if(W_DENS)
DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it);
- TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it);
- dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
- dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
- dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
+ dyni (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
+ Z_n_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_.*(KY==0)),Nx,Ny)));
+ Z_n_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_.*(KY==0)),Nx,Ny)));
+ end
+ if(W_TEMP)
+ TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it);
+ dyTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
+ dyTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
- end
+ Z_T_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny)));
+ Z_T_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_.*(KY==0)),Nx,Ny)));
+ end
end
end
-Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space
-
-for it = 1:numel(Ts5D)
- [~, it2D] = min(abs(Ts3D-Ts5D(it)));
- Np_i(:,:,it) = 0;
- for ij = 1:Nji
- Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
- Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it));
- end
-
- np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny)));
-end
% Post processing
disp('- post processing')
-% gyrocenter and particle flux from real space
-GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni drphi>
-GFlux_yi = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ni dzphi>
-GFlux_xe = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne drphi>
-GFlux_ye = zeros(1,Ns3D); % Gyrocenter flux Gamma = <ne dzphi>
-% Hermite energy spectrum
-epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D);
-epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D);
+% particle flux
+Gamma_x= zeros(Nx,Ny,Nz,Ns3D); % Radial particle transport
phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_maxx_avg = zeros(Nz,Ns3D); % Time evol. of the norm of phi
+phi_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
@@ -125,17 +119,16 @@ for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it))));
phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it))));
phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
+
+ if(W_DENS)
+ Gamma_x(:,:,iz,it) = dens_i(:,:,iz,it).*dyphi(:,:,iz,it);
+ end
- shear_maxx_maxy(iz,it) = max( max(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dr2phi(:,:,iz,it)))));
- shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dr2phi(:,:,iz,it)))));
+ shear_maxx_maxy(iz,it) = max( max(squeeze(-(dx2phi(:,:,iz,it)))));
+ shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dx2phi(:,:,iz,it)))));
+ shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dx2phi(:,:,iz,it)))));
+ shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dx2phi(:,:,iz,it)))));
- GFlux_xi(iz,it) = sum(sum(ni00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_yi(iz,it) = sum(sum(-ni00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_xe(iz,it) = sum(sum(ne00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly;
- GFlux_ye(iz,it) = sum(sum(-ne00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly;
-
Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn);
phi_kp_t(:,iz,it) = mean(Z_rth,2);
end
@@ -145,8 +138,6 @@ for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
[~, it2D] = min(abs(Ts3D-Ts5D(it)));
Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
- epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4);
- epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4);
end
%% Compute primary instability growth rate
diff --git a/matlab/setup.m b/matlab/setup.m
index 07937e16..b04e92c6 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -54,7 +54,8 @@ if (abs(CO) == 2) %Sugama operator
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
- INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5'''];
+% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5'''];
+ INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_10_J_5_N_20_kpm_8.0_NFLR_0.h5'''];
% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5'''];
elseif (CO == -1) % DGDK
disp('Warning, DGDK not debugged')
diff --git a/wk/HD_study.m b/wk/HD_study.m
index 7ec342c5..ca79206e 100644
--- a/wk/HD_study.m
+++ b/wk/HD_study.m
@@ -1,3 +1,5 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% DOUGHERTY
% bursts simulations
b_=[...
5e-1, 5e+0;...
@@ -7,8 +9,9 @@ b_=[...
1e-1, 1e-2;...
5e-2, 3e-3;...
5e-2, 2e-3;...
+ 1e-2, 1e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/
1e-2, 3e-3;...
- 1e-3, 2e-2;...
+ 1e-3, 2e-2;... % 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/
];
% converged turb plateau simulations
cp_=[...
@@ -16,8 +19,12 @@ cp_=[...
1e+0, 3e-3;...
1e+0, 2e-3;...
5e-1, 2e-3;...
+ 5e-1, 2e-2;...
1e-1, 2e-3;...
1e-1, 1e-3;...
+ 5e-2, 1e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02
+ 5e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-02_DGGK_mu_5e-04
+ 1e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
];
% moving/no turb plateau
dp_=[...
@@ -27,22 +34,21 @@ dp_=[...
1e-3, 2e-2;...
1e-3, 1e-2;...
1e-3, 5e-3;...
- 1e-3, 2.5e-3,...
+ 1e-3, 2.5e-3;...
+ 1e-3, 1e-3;... %v2.7_P_6_J_3/200x100_L_60_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_1e-03/
];
% not sure
rp_=[...
- 5e-2, 1e-3;...
- 1e-2, 1e-3;...
1e-2, 1e-4;...
- 1e-2, 5e-4;...
];
figure; set(gcf, 'Position', [100, 100, 900, 400])
+title('Hyperdiffusion study, Dougherty GK')
grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-xlabel('$\nu_{DGGK}$'); ylabel('$\mu_{HD}$'); hold on;
+xlabel('$\nu$'); ylabel('$\mu_{HD}$'); hold on;
% Trajectory of simulations
-
+if 0
% HD_study/200x100_L_200_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_1e-02/
mu_ = [0 1 2 5 5];
nu_ = [1 1 1 1 0.5];
@@ -54,9 +60,9 @@ nu_ = [1e-1 1e-1];
plot(nu_,mu_,'x--','DisplayName','N=300, L=100, P,J=2,1');
% HD_study/300x150_L_200_P_2_J_1_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_1e-02/
-% mu_ = [2e-3];
-% nu_ = [5e-1];
-% plot(nu_,mu_,'x--','DisplayName','N=300, L=100, P,J=2,1');
+mu_ = [2e-3 2e-2];
+nu_ = [5e-1 5e-1];
+plot(nu_,mu_,'x--','DisplayName','N=300, L=100, P,J=2,1');
% HD_study/100x50_L_50_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/
mu_ = [1e-2 5e-3 2e-3];
@@ -74,27 +80,109 @@ nu_ = [1e-2 1e-2 1e-2];
plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
% HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-02/
-mu_ = [3e-3 5e-4];
-nu_ = [1e-2 1e-2];
+mu_ = [3e-3 5e-4 1e-4];
+nu_ = [1e-2 1e-2 1e-2];
plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=4,2');
-
% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_5e-02_DGGK_CLOS_0_mu_3e-02/
mu_ = [3e-3 1e-3];
nu_ = [5e-2 5e-2];
plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
+% v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/
+mu_ = [2e-2];
+nu_ = [1e-1];
+% plot(nu_,mu_,'x--','DisplayName','N=200, L=120, P,J=6,3');
+end
scatter( b_(:,1), b_(:,2),'o',...
'MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
'DisplayName','Bursts');
scatter(cp_(:,1),cp_(:,2),'s',...
'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
'DisplayName','Converged Plateau');
-scatter(dp_(:,1),dp_(:,2),'v',...
+scatter(dp_(:,1),dp_(:,2),'d',...
'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
'DisplayName','Moving Plateau');
scatter(rp_(:,1),rp_(:,2),'h',...
'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
'DisplayName','not sure');
-% plot([1 1],[5e-5 5e-1],'--','Color',[0.4660 0.6740 0.1880]);
-legend('show','Location','Eastoutside')
\ No newline at end of file
+plot(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0], 'DisplayName','$\mu=0$');
+legend('show','Location','NorthWest')
+scatter(1,5e-5,80,'v','MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0]);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% SUGAMA
+% nu mu
+% bursts simulations
+b_=[...
+ 1e+0, 1.0e-2;... % 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/
+ 1e+0, 3.0e-3;... % 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/
+ 5e-1, 3.0e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02
+ 5e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02
+ 5e-1, 0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt
+ 1e-1, 2.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/
+ 1e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02/
+ 1e-2, 5.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
+ 1e-2, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
+ ];
+% converged turb plateau simulations
+cp_=[...
+ 1e-1, 0;... % v2.7_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_0e+00/
+ ];
+% moving/no turb plateau
+dp_=[...
+ 1e-1, 0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt
+ 1e-2, 1.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
+ 1e-2, 5.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
+ 1e-2, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
+ 1e-3, 1.0e-3;... % 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/
+ ];
+% not sure
+rp_=[...
+ 1e-1, 3.0e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02
+ 5e-2, 1.0e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02
+ ];
+figure; set(gcf, 'Position', [100, 100, 900, 400])
+title('Hyperdiffusion study, Sugama GK')
+grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('$\nu$'); ylabel('$\mu_{HD}$'); hold on;
+
+% Trajectory of simulations
+if 0
+% v2.7_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_0e+00/
+mu_ = [0.0];
+nu_ = [0.1];
+plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
+
+% Trajectory of simulations
+% v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/
+mu_ = [0.02];
+nu_ = [0.1];
+plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
+
+mu_ = [0.0];
+nu_ = [0.1];
+plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
+end
+
+scatter( b_(:,1), b_(:,2),'o',...
+ 'MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],'SizeData',80,...
+ 'DisplayName','Bursts');
+% scatter(cp_(:,1),cp_(:,2),'s',...
+% 'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
+% 'DisplayName','Converged Plateau');
+scatter(dp_(:,1),dp_(:,2),'d',...
+ 'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
+ 'DisplayName','Plateau');
+% scatter(rp_(:,1),rp_(:,2),'h',...
+% 'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
+% 'DisplayName','not sure');
+scatter(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],...
+ 'DisplayName','$\mu=0$');
+legend('show','Location','NorthWest')
+
+scatter(0.1,5e-5,80,'v','MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],...
+ 'DisplayName','$\mu=0$');
+scatter(0.5,5e-5,80,'v','MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],...
+ 'DisplayName','$\mu=0$');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 083fd6e3..63441138 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -6,10 +6,6 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
@@ -22,8 +18,8 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/test_3D_marconi/100x50x20_L_60_q0_1_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_2e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt';
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
@@ -41,14 +37,14 @@ default_plots_options
disp('Plots')
FMT = '.fig';
-if 1
+if 0
%% Time evolutions and growth rate
plot_time_evolution_and_gr
end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG = 5000; % Averaging time duration
+TAVG = 1000; % Averaging time duration
plot_radial_transport_and_shear
end
@@ -59,42 +55,50 @@ tstart = 0; tend = Ts3D(end); % time window
plot_space_time_diagrams
end
-if 0
-%% Averaged shear and Reynold stress profiles
-trange = its2D:ite2D;
-plot_shear_and_reynold_stress
-end
-
if 0
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
+% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
+tstart = 5000; tend = tstart+1000;
+% Chose the field to plot
+% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
+% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
+LOGSCALE = 1; TRENDS = 0;
plot_kperp_spectrum
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-t0 =00; iz = 1; ix = 1; iy = 1;
-skip_ = 1; DELAY = 1e-2*skip_;
+t0 =3000; iz = 1; ix = 1; iy = 1;
+skip_ =1; DELAY = 1e-2*skip_;
[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
-INTERP = 1; T = Ts3D; FRAMES = FRAMES_3D;
+INTERP = 0; T = Ts3D; FRAMES = FRAMES_3D;
% Field to plot
-% FIELD = dens_e; NAME = 'ne'; FIELDNAME = 'n_e';
-% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i';
-% FIELD = temp_e; NAME = 'Te'; FIELDNAME = 'n_i';
-% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
-% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
-FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-% Slice
+% FIELD = dens_e; NAME = 'ne'; FIELDNAME = 'n_e';
+% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i';
+% FIELD = dens_e-Z_n_e; NAME = 'ne_NZ';FIELDNAME = 'n_e^{NZ}';
+% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
+% FIELD = temp_e; NAME = 'Te'; FIELDNAME = 'n_i';
+% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
+% FIELD = temp_e-Z_T_e; NAME = 'Te_NZ';FIELDNAME = 'T_e^{NZ}';
+% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
+% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
+% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
+% FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
+
+% Sliced
% plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
% plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
-plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
+% plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
% Averaged
% plt = @(x) mean(x,1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
% plt = @(x) mean(x,2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
-% plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
+plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
FIELD = squeeze(plt(FIELD));
@@ -110,27 +114,32 @@ if 0
%% Photomaton : real space
% Chose the field to plot
-FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
-% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
-% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
-% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
-% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
-% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
-% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
+% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
+% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
+% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
+% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
+FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}';
+% FIELD = dens_i-Z_n_i; FNAME = 'ni_NZ'; FIELDLTX = 'n_i^{NZ}';
+% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
+% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
+% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
+% FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
+% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
+% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
% Chose when to plot it
-tf = [0 1 2 3];
+tf = 500:500:2500;
-% Slice
+% Sliced
ix = 1; iy = 1; iz = 1;
% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')']
% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')']
-% plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(z(iz))),')']
+plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(z=',num2str(round(z(iz)/pi)),'\pi)']
% Averaged
-% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_x']
-% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_y']
-plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_z']
+% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_x']
+% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_y']
+% plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_z']
%
@@ -150,5 +159,42 @@ plt_2 = @(x) x./max(max(x));
save_figure
end
+if 0
+%% Photomaton : k space
+
+% Chose the field to plot
+% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
+FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
+% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
+
+% Chose when to plot it
+tf = 500:500:2500;
+
+% Sliced
+ix = 1; iy = 1; iz = 1;
+% plt = @(x,it) abs(x(ix, :, :,it)); X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_x=',num2str(round(kx(ix))),')'];
+% plt = @(x,it) abs(x( :,iy, :,it)); X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_y=',num2str(round(ky(iy))),')'];
+plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)'];
+
+%
+TNAME = [];
+fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 300*numel(tf), 500])
+plt_2 = @(x) (fftshift(x,2));
+ for i_ = 1:numel(tf)
+ [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
+ subplot(1,numel(tf),i_)
+ DATA = plt_2(squeeze(plt(FIELD,it)));
+ pclr = pcolor(fftshift(X,2),fftshift(Y,2),DATA); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
+% colormap(bluewhitered); caxis([-1,1]);
+ xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
+ if(i_ > 1); set(gca,'ytick',[]); end;
+ title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
+ end
+ legend(latexize(FIELDLTX));
+save_figure
+end
+
%%
% ZF_fourier_analysis
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index 47f91c36..c3d7e300 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,9 +1,12 @@
%% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
%% Functions to modify preexisting fort.90 input file and launch on marconi
function [] = continue_run(outfilename)
- EXECNAME = 'helaz_2.8';
+ EXECNAME = 'helaz_3.1';
%% CLUSTER PARAMETERS
CLUSTER.PART = 'prod'; % dbg or prod
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
@@ -37,41 +40,41 @@ function [] = continue_run(outfilename)
fclose(fid);
% Find previous job2load
- if( numel(A{5}) == numel(' RESTART = .false.') )
+ if( numel(A{5}) == numel(' RESTART = .false.') )
A{5} = sprintf(' RESTART = .true.');
J2L = 0;
else
- line = A{35};
+ line = A{39};
line = line(end-2:end);
if(line(1) == '='); line = line(end); end;
J2L = str2num(line) + 1;
end
% Change job 2 load in fort.90
- A{35} = [' job2load = ',num2str(J2L)];
- disp(A{35})
+ A{39} = [' job2load = ',num2str(3)];
+ disp(A{39})
% Change time step
- A{3} = [' dt = 0.005'];
+ A{3} = [' dt = 0.01'];
% Increase endtime
A{4} = [' tmax = 10000'];
% Put non linear term back
- A{41} = [' NL_CLOS = -1'];
+ A{45} = [' NL_CLOS = -1'];
% change HD
- line_= A{43};
+ line_= A{47};
mu_old = str2num(line_(13:end));
- A{43} = [' mu = ',num2str(mu_old*2)];
+ A{47} = [' mu = ',num2str(mu_old)];
% change L
line_= A{14};
- L_old = str2num(line_(8:end));
- A{14} = [' Lx = ',num2str(L_old)];
- A{16} = [' Ly = ',num2str(L_old)];
+ L_old = str2num(line_(12:end));
+ A{14} = [' Lx = ',num2str(L_old)];
+ A{16} = [' Ly = ',num2str(L_old)];
% change eta N
- line_= A{53};
+ line_= A{57};
etan_old = str2num(line_(13:end));
- A{53} = [' eta_n = ',num2str(etan_old)];
+ A{57} = [' eta_n = ',num2str(etan_old)];
% change eta B
- line_= A{55};
+ line_= A{59};
etab_old = str2num(line_(13:end));
- A{55} = [' eta_B = ',num2str(etab_old)];
+ A{59} = [' eta_B = ',num2str(etab_old)];
% 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 c0ac260d..efbbcce1 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -25,7 +25,7 @@ TMAX = 200; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/50; % Sampling per time unit for 5D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 00;
@@ -141,8 +141,9 @@ end
if 1
%% Parameter scan over CO
-PMAXE = 6; PMAXI = 6;
-JMAXE = 3; JMAXI = 3;
+P=2; J=1;
+N = 20; % Frequency gridpoints (Nkx = N/2)
+L = 75; % Size of the squared frequency domain
TMAX = 200;
DT = 0.01;
CO_A = [4];
@@ -151,7 +152,8 @@ Nparam = numel(CO_A);
param_name = 'CO';
ETAN = 2.0;
NU = 1e-1; % Collision frequency
-
+PMAXE = P; PMAXI = P;
+JMAXE = J; JMAXI = J;
Bohm_transport = zeros(Nparam,1);
gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
gamma_Nipj = zeros(Nparam,floor(N/2)+1);
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 85f76881..31afde74 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,3 +1,6 @@
addpath(genpath('../matlab')) % ... add
%% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/test_3D_marconi/100x50_L_60_P_2_J_1_eta_Inf_nu_1e-01_DGGK_mu_0e+00/out.txt');
\ No newline at end of file
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index f6f67373..5bb94669 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -8,21 +8,21 @@ NU = 0.1; % Collision frequency
ETAN = 1.0/0.6; % Density gradient drive (R/Ln)
NU_HYP = 1.0;
%% GRID AND GEOMETRY PARAMETERS
-N = 50; % Frequency gridpoints (Nkx = N/2)
-L = 30; % Size of the squared frequency domain
-Nz = 10; % number of perpendicular planes (parallel grid)
+N = 150; % Frequency gridpoints (Nkx = N/2)
+L = 100; % Size of the squared frequency domain
+Nz = 1; % number of perpendicular planes (parallel grid)
q0 = 1.0; % safety factor
shear = 0.0; % magnetic shear
eps = 0.0; % inverse aspect ratio
-P = 2;
-J = 1;
+P = 4;
+J = 2;
%% TIME PARAMETERS
-TMAX = 500; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 1000; % Maximal time unit
+DT = 1e-3; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 2; % Sampling per time unit for 3D arrays
-SPS5D = 1; % Sampling per time unit for 5D arrays
+SPS3D = 1/2; % Sampling per time unit for 3D arrays
+SPS5D = 1/20; % 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;
@@ -32,10 +32,10 @@ JOB2LOAD= 0;
CO = 1;
CLOS = 0; % Closure model (0: =0 truncation)
NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID = 'HD_study'; % Name of the simulation
-SIMID = 'test_3D'; % Name of the simulation
+SIMID = 'HD_study'; % Name of the simulation
+% SIMID = 'test_3D'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
-NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
+NON_LIN = -1; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
INIT_ZF = 0; ZF_AMP = 0.0;
INIT_BLOB = 0; WIPE_TURB = 0;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index f93e63cf..1c4f2ac0 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -18,36 +18,36 @@ NP_P = 2; % MPI processes along p
NP_KX = 24; % MPI processes along kx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
ETAN = 1.0/0.6; % Density gradient drive (R/Ln)
NU_HYP = 1.0;
%% GRID PARAMETERS
-N = 100; % Frequency gridpoints (Nkx = N/2)
-L = 60; % Size of the squared frequency domain
-Nz = 20; % number of perpendicular planes (parallel grid)
+N = 150; % Frequency gridpoints (Nkx = N/2)
+L = 100; % Size of the squared frequency domain
+Nz = 1; % number of perpendicular planes (parallel grid)
q0 = 1.0; % q factor ()
shear = 0.0; % magnetic shear
eps = 0.0; % inverse aspect ratio
-P = 2;
-J = 1;
+P = 4;
+J = 2;
%% TIME PARAMETERS
-TMAX = 200; % Maximal time unit
+TMAX = 10000; % 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
SPS3D = 1/2; % Sampling per time unit for 3D arrays
-SPS5D = 1/100; % Sampling per time unit for 5D arrays
+SPS5D = 1/500; % 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;
+RESTART = 1; % To restart from last checkpoint
+JOB2LOAD= 3;
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO = 1;
+CO = 2;
CLOS = 0; % Closure model (0: =0 truncation)
-NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'test_3D_marconi'; % Name of the simulation
-% SIMID = 'HD_study'; % Name of the simulation
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+% SIMID = 'test_3D_marconi'; % Name of the simulation
+SIMID = 'HD_study'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index b81ce091..d320b533 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -65,14 +65,16 @@ subplot(224)
%% FCGK
P_ = 6; J_ = 3;
-mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5';
-kp = 1.2; matidx = round(kp/(8/150));
+mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5';
+kp = 1.0; matidx = round(kp/(8/150));
+% matidx = sprintf('%5.5i',matidx);
+% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_1_kpm_4.0_NFLR_20_m_0.h5';
matidx = sprintf('%5.5i',matidx);
FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
FCGK_ie = h5read(mat_file_name,['/',matidx,'/Ciepj/CiepjF'])+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
figure
-MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=20,J=10, k=',num2str(kp)]);
+MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=6,J=3, k=',num2str(kp),' (',matidx,')']);
% MAT = 1i*FCGK_ei; suptitle('FCGK ei,P=20,J=10');
% MAT = 1i*FCGK_ie; suptitle('FCGK ie,P=20,J=10');
subplot(221)
@@ -88,39 +90,56 @@ subplot(224)
imagesc(imag(MAT)<0);
title('imag$<$0');
-%% Eigenvalue analysis
+%% Single eigenvalue analysis
+
+% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_2_J_1_N_1_kpm_4.0_NFLR_5.h5';
+% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb.NFLR/k.4/self.NFLR.8.0.h5';
+matidx = 74;
+
+matidx = sprintf('%5.5i',matidx);disp(matidx);
+
+% MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
+MAT = h5read(mat_file_name,['/Caapj/Ciipj']);
+
+gmax = max(real(eig(MAT)));
+
+wmax = max(imag(eig(MAT)));
+figure
+imagesc((MAT)); colormap(bluewhitered)
+title(['$\gamma_{max}=',num2str(gmax),'$'])
+%% Eigenvalue spectrum analysis
if 0
%%
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_10_J_5_N_20_kpm_8.0_NFLR_0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_2_J_1_N_20_kpm_8.0_NFLR_5.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5',...
};
-CONAME_A = {'SG','PA','FC NFLR 4', 'FC NFLR k2'};
-kp_a = linspace(0,8,149);
-gmax = zeros(size(kp_a));
-wmax = zeros(size(kp_a));
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 10 5 NFLR 0', 'FC 2 1 NFLR 5', 'FC 6 3 NFLR 0', 'FC 6 3 NFLR 8'};
figure
-for j_ = 1:4
- mat_file_name = mfns{j_};
- i = 1;
- for kp = kp_a
- matidx = floor(kp/(8/149));
- matidx = sprintf('%5.5i',matidx);disp(matidx)
+for j_ = 1:numel(mfns)
+ mat_file_name = mfns{j_}; disp(mat_file_name);
+ kp_a = h5read(mat_file_name,'/coordkperp');
+ gmax = zeros(size(kp_a));
+ wmax = zeros(size(kp_a));
+ for idx_ = 0:numel(kp_a)-1
+ matidx = sprintf('%5.5i',idx_);
MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
- gmax(i) = max(real(eig(MAT)));
- wmax(i) = max(imag(eig(MAT)));
- i = i + 1;
+ gmax(idx_+1) = max(real(eig(MAT)));
+ wmax(idx_+1) = max(imag(eig(MAT)));
end
- subplot(121)
+% subplot(121)
plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
- subplot(122)
- plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
+% subplot(122)
+% plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
end
- subplot(121)
+% subplot(121)
legend('show'); grid on;
xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)')
- subplot(122)
-legend('show'); grid on;
-xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)')
+% subplot(122)
+% legend('show'); grid on;
+% xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)')
end
--
GitLab