diff --git a/matlab/compile_results.m b/matlab/compile_results.m index b79637a..44177cf 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,139 +1,139 @@ CONTINUE = 1; 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 Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; GGAMMA_ = []; PGAMMA_ = []; PHI_ = []; DENS_E_ = []; DENS_I_ = []; TEMP_E_ = []; TEMP_I_ = []; Ts0D_ = []; Ts3D_ = []; 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) % Load results of simulation #JOBNUM load_results % 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 tmp = Nipj_; sz = size(tmp); Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; tmp = Nepj_; sz = size(tmp); Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; end if W_SAPJ tmp = Sipj_; sz = size(tmp); Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; tmp = Sepj_; sz = size(tmp); 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 GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); Ts0D_ = cat(1,Ts0D_,Ts0D); end if W_PHI || W_NA00 Ts3D_ = cat(1,Ts3D_,Ts3D); end if W_PHI PHI_ = cat(4,PHI_,PHI); end if W_NA00 Ni00_ = cat(4,Ni00_,Ni00); Ne00_ = cat(4,Ne00_,Ne00); end if W_DENS DENS_E_ = cat(4,DENS_E_,DENS_E); DENS_I_ = cat(4,DENS_I_,DENS_I); end if W_TEMP TEMP_E_ = cat(4,TEMP_E_,TEMP_E); TEMP_I_ = cat(4,TEMP_I_,TEMP_I); end if W_NAPJ || W_SAPJ Ts5D_ = cat(1,Ts5D_,Ts5D); end if W_NAPJ Nipj_ = cat(6,Nipj_,Nipj); Nepj_ = cat(6,Nepj_,Nepj); end if W_SAPJ Sipj_ = cat(6,Sipj_,Sipj); Sepj_ = cat(6,Sepj_,Sepj); end % Evolution of simulation parameters load_params 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]; JOBFOUND = JOBFOUND + 1; LASTJOB = JOBNUM; Pe_old = Pe_new; Je_old = Je_new; Pi_old = Pi_new; Ji_old = Ji_new; elseif (JOBNUM > JOBNUMMAX) CONTINUE = 0; disp(['found ',num2str(JOBFOUND),' results']); end JOBNUM = JOBNUM + 1; end GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_; Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_; Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts3D = Ts3D_; 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_; 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 d333c1b..1253d90 100644 --- a/matlab/dnjs.m +++ b/matlab/dnjs.m @@ -1,28 +1,28 @@ 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); Factar(1) = 1.0; f_ = 1; for ii_ = 2:(n+j+s)+1 f_ = (ii_-1) * f_; Factar(ii_) = f_; end result = 0; for il3 = 1:numel(L3) for il2 = 1:numel(L2) for il1 = 1:numel(L1) result = result + Factar(il1 + il2 + il3 -2)... *L1(il1) * L2(il2) * L3(il3); end end end end \ No newline at end of file diff --git a/matlab/dnjs_fact.m b/matlab/dnjs_fact.m new file mode 100644 index 0000000..ee5faf8 --- /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 0000000..013ac5f --- /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 f7eede8..400e63c 100644 --- a/matlab/plot_kernels.m +++ b/matlab/plot_kernels.m @@ -1,46 +1,216 @@ +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; nmax2=10; kmax=7; 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 2370dd0..238faa0 100644 --- a/matlab/plots/plot_kperp_spectrum.m +++ b/matlab/plots/plot_kperp_spectrum.m @@ -1,37 +1,36 @@ [~,itstart] = min(abs(Ts3D-tstart)); [~,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)); kp_ip = kx; [thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip); [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)),' -GFlux_yi = zeros(1,Ns3D); % Gyrocenter flux Gamma = -GFlux_xe = zeros(1,Ns3D); % Gyrocenter flux Gamma = -GFlux_ye = zeros(1,Ns3D); % Gyrocenter flux Gamma = -% 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 shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % . % Kperp spectrum interpolation %full kperp points kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]); % interpolated kperps nk_noAA = floor(2/3*numel(kx)); kp_ip = kx; [thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip); [xn,yn] = pol2cart(thg,rg); [ky_s, sortIdx] = sort(ky); [xc,yc] = meshgrid(ky_s,kx); phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D); % for it = 1:numel(Ts3D) % Loop over 2D aX_XYays for iz = 1:numel(z) NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it)))); 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 end % 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 disp('- growth rate') % Find max value of transport (end of linear mode) [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2); [~,itmax] = min(abs(Ts3D-tmax)); tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax); [~,its3D_lin] = min(abs(Ts3D-tstart)); [~,ite3D_lin] = min(abs(Ts3D-tend)); g_I = zeros(Nkx,Nky,Nz); for ikx = 1:Nkx for iky = 1:Nky for iz = 1:Nz [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); end end end [gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3); kmax_I = abs(ky(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/Nx/Ny)^2); % [~,itmax] = min(abs(Ts2D-tmax)); % tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax); [~,its3D_lin] = min(abs(Ts3D-tstart)); [~,ite3D_lin] = min(abs(Ts3D-tend)); g_II = zeros(Nkx,Nky); for ikx = 1:Nkx for iky = 1 for iz = 1:Nz [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); end end end [gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3); kmax_II = abs(kx(ikmax_II)); diff --git a/matlab/setup.m b/matlab/setup.m index 07937e1..b04e92c 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,150 +1,151 @@ %% ________________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nx = N; % x grid resolution GRID.Lx = L; % x length GRID.Ny = N * (1-KXEQ0) + KXEQ0; % y '' GRID.Ly = L * (1-KXEQ0); % y '' GRID.Nz = Nz; % z resolution GRID.q0 = q0; % q factor GRID.shear = shear; % shear GRID.eps = eps; % inverse aspect ratio % Model parameters MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end; MODEL.mu = MU; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % mass ratio sqrt(m_a/m_i) MODEL.sigma_e = 0.0233380; MODEL.sigma_i = 1.0; % charge q_a/e MODEL.q_e =-1.0; MODEL.q_i = 1.0; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.eta_n = ETAN; % source term kappa for HW MODEL.eta_T = ETAT; % Temperature MODEL.eta_B = ETAB; % Magnetic MODEL.lambdaD = LAMBDAD; % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KX0KH),'_A_',num2str(A0KH)]; end; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end; INITIAL.INIT_ZF = INIT_ZF; if (WIPE_TURB); INITIAL.wipe_turb = '.true.'; else; INITIAL.wipe_turb = '.false.';end; if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end; INITIAL.init_background = (INIT_ZF>0)*ZF_AMP; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; INITIAL.mat_file = '''null'''; if (abs(CO) == 2) %Sugama operator INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''']; 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') end % Naming and creating input file switch abs(CO) case 0; CONAME = 'LB'; case 1; CONAME = 'DG'; case 2; CONAME = 'SG'; case 3; CONAME = 'PA'; case 4; CONAME = 'FC'; otherwise; CONAME ='UK'; end if (CO <= 0); CONAME = [CONAME,'DK']; else; CONAME = [CONAME,'GK']; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end if (PMAXE == PMAXI) && (JMAXE == JMAXI) degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; else degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; end degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',... CONAME,'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end if (Nz == 1) resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'_']; gridname = ['L_',num2str(L),'_']; else resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'x',num2str(GRID.Nz),'_']; gridname = ['L_',num2str(L),'_q0_',num2str(q0),'_']; end 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.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end; OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT); if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end; if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.false.';end; OUTPUTS.resfile0 = '''outputs'''; OUTPUTS.rstfile0 = '''checkpoint'''; OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) 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; MAKE = 'cd ..; make; cd wk'; system(MAKE); %% disp(['Set up ',SIMID]); disp([resolution,gridname,degngrad]); if RESTART disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/wk/HD_study.m b/wk/HD_study.m index 7ec342c..ca79206 100644 --- a/wk/HD_study.m +++ b/wk/HD_study.m @@ -1,100 +1,188 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% DOUGHERTY % bursts simulations b_=[... 5e-1, 5e+0;... 1e-1, 1e-1;... 1e-1, 3e-2;... 1e-1, 2e-2;... 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_=[... 1e+0, 1e-2;... 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_=[... 1e+0, 5e+0;... 1e+0, 2e+0;... 1e+0, 1e+0;... 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]; plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1'); % HD_study/300x150_L_200_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/ mu_ = [1e-2 1e-3]; 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]; nu_ = [1e-1 1e-1 1e-1]; plot(nu_,mu_,'x-','DisplayName','N=100, L=050, P,J=2,1'); % HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02/ mu_ = [3e-2 1e-3 0]; nu_ = [1e-1 1e-1 1e-1]; plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1'); % HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-02/ mu_ = [3e-3 1e-3 1e-4]; 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 083fd6e..6344113 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,154 +1,200 @@ addpath(genpath('../matlab')) % ... add addpath(genpath('../matlab/plots')) % ... add %% Directory of the simulation if 1% Local results 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,'/']; CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD); system(CMD); else% Marconi results outfile =''; outfile =''; 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 %% Load the results % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 20; compile_results %Compile the results from first output found to JOBNUMMAX if existing %% Post-processing post_processing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 if 0 %% Space time diagramms cmax = 0.01 % max of the colorbar for transport 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)); % Naming GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS]; % Create movie (gif or mp4) create_gif % create_mov end 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'] % TNAME = []; fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350]) plt_2 = @(x) x./max(max(x)); 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((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap(bluewhitered); caxis([-1,1]); xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); end legend(latexize(FIELDLTX)); 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 47f91c3..c3d7e30 100644 --- a/wk/continue_multiple_runs_marconi.m +++ b/wk/continue_multiple_runs_marconi.m @@ -1,89 +1,92 @@ %% 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 if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end; CLUSTER.MEM = '64GB'; % Memory CLUSTER.JNAME = 'HeLaZ';% Job name NP_P = 2; % MPI processes along p NP_KX = 24; % MPI processes along kx % Compute processes distribution Ntot = NP_P * NP_KX; Nnodes = ceil(Ntot/48); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kx CLUSTER.CPUPT = '1'; % CPU per task %% RESDIR = ['../',outfilename(46:end-8),'/']; BASIC.RESDIR = RESDIR; FORT90 = [RESDIR,'fort.90']; % Read txt into cell A fid = fopen(FORT90,'r'); i = 1; tline = fgetl(fid); A{i} = tline; while ischar(tline) i = i+1; tline = fgetl(fid); A{i} = tline; end 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) if A{i+1} == -1 fprintf(fid,'%s', A{i}); break else fprintf(fid,'%s\n', A{i}); end end % Copy fort.90 into marconi write_sbash_marconi % Launch the job system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); end \ No newline at end of file diff --git a/wk/linear_study.m b/wk/linear_study.m index c0ac260..efbbcce 100644 --- a/wk/linear_study.m +++ b/wk/linear_study.m @@ -1,224 +1,226 @@ %clear all; addpath(genpath('../matlab')) % ... add default_plots_options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS % NU = 1.0; % Collision frequency TAU = 1.0; % e/i temperature ratio ETAB = 1.0; % ETAN = ETAN; % Density gradient ETAT = 0.0; % Temperature gradient NU_HYP = 0.0; % Hyperdiffusivity coefficient LAMBDAD = 0.0; NOISE0 = 1.0e-5; %% GRID PARAMETERS N = 50; % Frequency gridpoints (Nkx = N/2) L = 75; % Size of the squared frequency domain KXEQ0 = 1; % put kx = 0 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 PARMETERS 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; %% OPTIONS SIMID = 'v3.1_lin_analysis'; % Name of the simulation NON_LIN = 0 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : 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) INIT_PHI= 0; % Start simulation with a noisy phi INIT_ZF = 0; % Start simulation with a noisy phi %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 0; W_PHI = 0; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; W_DENS = 0; W_TEMP = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused % DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1) JOBNUM = 00; KPAR = 0.0; % Parellel wave vector component HD_CO = 0.5; % Hyper diffusivity cutoff ratio kmax = N*pi/L;% Highest fourier mode MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient Nz = 1; % number of perpendicular planes (parallel grid) q0 = 1.0; % safety factor shear = 0.0; % magnetic shear eps = 0.0; % inverse aspect ratio %% PARAMETER SCANS if 0 %% Parameter scan over PJ % PA = [2 4 6 10]; % JA = [1 2 3 5]; PA = [6]; JA = [3]; DTA= DT*ones(size(JA));%./sqrt(JA); % DTA= DT; mup_ = MU_P; muj_ = MU_J; Nparam = numel(PA); 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,floor(SPS2D*TMAX)); for i = 1:Nparam % Change scan parameter PMAXE = PA(i); PMAXI = PA(i); JMAXE = JA(i); JMAXI = JA(i); DT = DTA(i); setup % Run linear simulation % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz_3 1 1; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 2 3; cd ../../../wk']) % Load and process results %% filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results tend = Ts3D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); for ikx = 1:N/2+1 gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend)))))); Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:))); end tend = Ts5D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts5D-tstart)); [~,itend] = min(abs(Ts5D-tend)); for ikx = 1:N/2+1 gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2))); end gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0)); gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0)); % kymax = abs(kx(ikymax)); % Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; % Clean output system(['rm -r ',BASIC.RESDIR]); end if 1 %% Plot SCALE = 1;%sqrt(2); fig = figure; FIGNAME = 'linear_study'; plt = @(x) x; % subplot(211) 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); semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),... 'Color',clr,... 'LineStyle',linestyle{1},'Marker','^',... 'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\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(kx)]); title(['$\eta=',num2str(ETAB/ETAN),'$, $\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 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]; CONAME_A = {}; 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); for i = 1:Nparam % Change scan parameter CO = CO_A(i); setup CONAME_A{i} = CONAME; % Run linear simulation system(... ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']... ) %% Load an process results filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results tend = Ts3D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); for ikx = 1:N/2+1 gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend)))))); Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:))); end tend = Ts5D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts5D-tstart)); [~,itend] = min(abs(Ts5D-tend)); for ikx = 1:N/2+1 gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2))); end gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0)); gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0)); % kymax = abs(kx(ikymax)); % Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; % Clean output system(['rm -r ',BASIC.RESDIR]); end if 1 %% Plot SCALE = 1;%sqrt(2); fig = figure; FIGNAME = 'linear_study'; plt = @(x) x; % subplot(211) 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); semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),... 'Color',clr,... 'LineStyle',linestyle{1},'Marker','^',... 'DisplayName',[CONAME_A{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(kx)]); title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu=',num2str(NU),'$',', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) 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 if 0 %% Plot fig = figure; FIGNAME = 'mixing_length'; plot(eta_B, Bohm_transport) grid on; xlabel('$L_n/L_B$'); ylabel('$\eta\gamma_{max}/k_{max}^2$'); title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']); end %% end \ No newline at end of file diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m index 85f7688..31afde7 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 f6f6737..5bb9466 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,78 +1,78 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS 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; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK) 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; %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; W_DENS = 1; W_TEMP = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio % kmaxcut = 2.5; MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB) INIT_PHI= 1; % Start simulation with a noisy phi and moments MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% Setup and file management setup system('rm fort.90'); outfile = [BASIC.RESDIR,'out.txt']; disp(outfile); diff --git a/wk/marconi_run.m b/wk/marconi_run.m index f93e63c..1c4f2ac 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,106 +1,106 @@ clear all; addpath(genpath('../matlab')) % ... add SUBMIT = 1; % To submit the job automatically % EXECNAME = 'helaz_dbg'; EXECNAME = 'helaz_3.1'; for ETAN = [1/0.6] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS 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 CLUSTER.JNAME = 'HeLaZ';% Job name 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 INIT_ZF = 0; ZF_AMP = 0.0; INIT_BLOB = 0; WIPE_TURB = 0; %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; W_DENS = 1; W_TEMP = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio % kmaxcut = 2.5; MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB) INIT_PHI= 1; % Start simulation with a noisy phi and moments MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f % Compute processes distribution Ntot = NP_P * NP_KX; Nnodes = ceil(Ntot/48); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kx CLUSTER.CPUPT = '1'; % CPU per task %% Run file management scripts setup write_sbash_marconi system('rm fort.90 setup_and_run.sh batch_script.sh'); if(mod(NP_P*NP_KX,48)~= 0) disp('WARNING : unused cores (ntot cores must be a 48 multiple)'); end if(SUBMIT) system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); end disp('done'); end \ No newline at end of file diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m index b81ce09..d320b53 100644 --- a/wk/plot_cosol_mat.m +++ b/wk/plot_cosol_mat.m @@ -1,126 +1,145 @@ % MAT = results.iCa; % figure % suptitle('FCGK,P=6,J=3'); % subplot(221) % imagesc(log(abs(MAT))); % title('log abs') % subplot(222) % imagesc(imag(MAT)); colormap(bluewhitered) % title('imag'); colorbar % subplot(223) % imagesc(imag(MAT)>0); % title('imag$>$0'); % subplot(224) % imagesc(imag(MAT)<0); % title('imag$<$0'); % % %% SGGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; SGGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj'); SGGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); SGGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=0'); % MAT = 1i*SGGK_ei; suptitle('SGGK ei,P=20,J=10'); % MAT = 1i*SGGK_ie; suptitle('SGGK ie,P=20,J=10'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% PAGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; PAGK_self = h5read(mat_file_name,'/00000/Caapj/Ceepj'); % PAGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); % PAGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=0'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% 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) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); 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