diff --git a/.gitignore b/.gitignore index 491f19f..0b32f91 100644 --- a/.gitignore +++ b/.gitignore @@ -1,32 +1,33 @@ .vimrc *.swp *~ *.asv *.m.orig *.DS_Store *.fig *.pdf *.eps *.mov *.mp4 *.emf *.pptx *.jpg *.jpeg *.mat *.dat *.mod *.h5 *.o *.a +*.png logs/ results/ results_old/ mod/ obj/ bin/ .gitignore wk/fort.90 .directory checkpoint/ FM/ diff --git a/matlab/create_gif.m b/matlab/create_gif.m index 826ea6d..f2445a7 100644 --- a/matlab/create_gif.m +++ b/matlab/create_gif.m @@ -1,56 +1,56 @@ title1 = GIFNAME; FIGDIR = BASIC.RESDIR; GIFNAME = [FIGDIR, GIFNAME,'.gif']; % Set colormap boundaries hmax = max(max(max(FIELD))); hmin = min(min(min(FIELD))); flag = 0; if hmax == hmin disp('Warning : h = hmin = hmax = const') else % Setup figure frame fig = figure('Color','white','Position', [100, 100, 400, 400]); pcolor(X,Y,FIELD(:,:,1)); % to set up - colormap jet + colormap gray axis tight manual % this ensures that getframe() returns a consistent size if INTERP shading interp; end in = 1; nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); for n = FRAMES % loop over selected frames scale = max(max(abs(FIELD(:,:,n)))); pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot if INTERP shading interp; end set(pclr, 'edgecolor','none'); axis square; % caxis([min(min(FIELD(:,:,n))),max(max(FIELD(:,:,n)))]); xlabel(XNAME); ylabel(YNAME); %colorbar; title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... ,', scaling = ',sprintf('%.1e',scale)]); drawnow % Capture the plot as an image frame = getframe(fig); im = frame2im(frame); [imind,cm] = rgb2ind(im,32); % Write to the GIF File if in == 1 imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf); else imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY); end % terminal info while nbytes > 0 fprintf('\b') nbytes = nbytes - 1; end nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:))); in = in + 1; end disp(' ') disp(['Gif saved @ : ',GIFNAME]) end diff --git a/matlab/create_gif_3D.m b/matlab/create_gif_3D.m index 188349d..646080f 100644 --- a/matlab/create_gif_3D.m +++ b/matlab/create_gif_3D.m @@ -1,49 +1,48 @@ title1 = GIFNAME; FIGDIR = BASIC.RESDIR; GIFNAME = [FIGDIR, GIFNAME,'.gif']; % Setup figure frame fig = figure('Color','white','Position', [100, 100, 400, 400]); plt = @(x) squeeze(reshape(x(:,:,1),[],1)); scale_x = max(abs(plt(X(:,:,1)))); scale_y = max(abs(plt(Y(:,:,1)))); scale_z = max(abs(plt(Z(:,:,1)))); plot3(plt(X)/scale_x,plt(Y)/scale_y,plt(Z)/scale_z,'.k','MarkerSize',MARKERSIZE); view(VIEW); - colormap jet axis tight manual % this ensures that getframe() returns a consistent size in = 1; nbytes = fprintf(2,'frame %d/%d',in,numel(FRAMES)); for n = FRAMES % loop over selected frames plt = @(x) squeeze(reshape(x(:,:,n),[],1)); scale_x = max(abs(plt(X))); scale_y = max(abs(plt(Y))); scale_z = max(abs(plt(Z))); plot3(plt(X)/scale_x,plt(Y)/scale_y,plt(Z)/scale_z,'.k','MarkerSize',MARKERSIZE); view(VIEW); xlabel(XNAME); ylabel(YNAME); zlabel(ZNAME); title(['$t \approx$', sprintf('%.3d',ceil(T(n)))... ,', scaling = ',sprintf('%.1e',scale_x)]); grid on; drawnow % Capture the plot as an image frame = getframe(fig); im = frame2im(frame); [imind,cm] = rgb2ind(im,32); % Write to the GIF File if in == 1 imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf); else imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY); end % terminal info while nbytes > 0 fprintf('\b') nbytes = nbytes - 1; end nbytes = fprintf(2,'frame %d/%d',n,numel(FRAMES)); in = in + 1; end disp(' ') disp(['Gif saved @ : ',GIFNAME]) diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m index ed53cc0..56e500d 100644 --- a/matlab/default_plots_options.m +++ b/matlab/default_plots_options.m @@ -1,9 +1,10 @@ %% Default values for plots set(0,'defaulttextInterpreter','latex') set(0, 'defaultAxesTickLabelInterpreter','latex'); set(0, 'defaultLegendInterpreter','latex'); set(0,'defaultAxesFontSize',16) set(0, 'DefaultLineLineWidth', 2.0); line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];... [0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]]; -line_styles = {'-','-.','--',':'}; \ No newline at end of file +line_styles = {'-','-.','--',':'}; +marker_styles = {'^','v','o','s'}; \ No newline at end of file diff --git a/matlab/save_figure.m b/matlab/save_figure.m index e71ee64..c5a9c72 100644 --- a/matlab/save_figure.m +++ b/matlab/save_figure.m @@ -1,5 +1,8 @@ %% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME) % and parameters -FIGNAME = [BASIC.RESDIR, FIGNAME,FMT]; -saveas(fig,FIGNAME); -disp(['Figure saved @ : ',FIGNAME]) \ No newline at end of file +if ~exist([BASIC.RESDIR,'/fig'], 'dir') + mkdir([BASIC.RESDIR,'/fig']) +end +saveas(fig,[BASIC.RESDIR,'/fig/', FIGNAME,'.fig']); +saveas(fig,[BASIC.RESDIR, FIGNAME,'.png']); +disp(['Figure saved @ : ',[BASIC.RESDIR, FIGNAME,'.png']]) \ No newline at end of file diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index 1ffdb66..6c739bb 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,344 +1,345 @@ %% Load results +% JOBNUM = 0; load_results; compile_results %% 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); Ns5D = numel(Ts5D); Ns2D = numel(Ts2D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts2D = Ts2D'; if strcmp(OUTPUTS.write_non_lin,'.true.') Si00 = squeeze(Sipj(1,1,:,:,:)); Se00 = squeeze(Sepj(1,1,:,:,:)); end %% Build grids Nkr = numel(kr); Nkz = numel(kz); [KZ,KR] = meshgrid(kz,kr); Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz); dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1); Lk = max(Lkr,Lkz); dr = 2*pi/Lk; dz = 2*pi/Lk; Nr = max(Nkr,Nkz); Nz = Nr; r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r); z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z); [ZZ,RR] = meshgrid(z,r); %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nr,Nz,Ns2D); ni00 = zeros(Nr,Nz,Ns2D); si00 = zeros(Nr,Nz,Ns5D); phi = zeros(Nr,Nz,Ns2D); drphi = zeros(Nr,Nz,Ns2D); dzphi = zeros(Nr,Nz,Ns2D); for it = 1:numel(Ts2D) NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz))); ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz))); phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz))); drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz))); dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz))); end if strcmp(OUTPUTS.write_non_lin,'.true.') for it = 1:numel(Ts5D) si00(:,:,it) = real(fftshift(ifft2(squeeze(Si00(:,:,it)),Nr,Nz))); end end % Post processing disp('- post processing') E_pot = zeros(1,Ns2D); % Potential energy n^2 E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2 ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi| Flux_ri = zeros(1,Ns2D); % Particle flux Gamma = Flux_zi = zeros(1,Ns2D); % Particle flux Gamma = Flux_re = zeros(1,Ns2D); % Particle flux Gamma = Flux_ze = zeros(1,Ns2D); % Particle flux Gamma = Ne_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj Ni_norm = zeros(Npi,Nji,Ns5D);% . if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj Si_norm = zeros(Npi,Nji,Ns5D);% . Sne00_norm = zeros(1,Ns2D); % Time evol. of the amp of e nonlin term Sni00_norm = zeros(1,Ns2D); % end Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2; for it = 1:numel(Ts2D) % Loop over 2D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); E_pot(it) = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id E_kin(it) = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr; ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz)))); Flux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; Flux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; end E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); dEdt = diff(E_pot+E_kin)./dt2D; for it = 1:numel(Ts5D) % Loop over 5D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz; if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz; Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz; Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz; end end %% Compute growth rate disp('- growth rate') tend = Ts2D(end); tstart = 0.6*tend; g_ = zeros(Nkr,Nkz); [~,ikr0KH] = min(abs(kr-KR0KH)); for ikr = 1:Nkr for ikz = 1:Nkz g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend); end end % gkr0kz_Ni00 = max(real(g_(:,:)),[],1); gkr0kz_Ni00 = real(g_(ikr0KH,:)); %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Plots') FMT = '.fig'; %% Time evolutions fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)]; subplot(221); for ip = 1:Npe for ij = 1:Nje plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$']; semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on; end end grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); subplot(222) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$']; semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on; end end grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); subplot(223) plot(Ts2D,Flux_ri,'-','DisplayName','$\Gamma_{ri}$'); hold on; plot(Ts2D,Flux_zi,'-','DisplayName','$\Gamma_{zi}$'); hold on; plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$') plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$') grid on; xlabel('$t$'); ylabel('$\Gamma$'); %legend('show'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$S_i^{',num2str(ip-1),num2str(ij-1),'}$']; semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on; end end grid on; xlabel('$t$'); ylabel('$S$'); %legend('show'); else %% Growth rate subplot(224) [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz(1:ikr0KH)/kr(ikr0KH),gkr0kz_Ni00(1:ikr0KH)/(KR0KH*A0KH),... 'DisplayName',['J = ',num2str(JMAXI)]); xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$'); xlim([0. 1.0]); %ylim([0. 0.04]); end save_figure %% if 1 %% Show frame in real space -tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf)); +tf = 1000; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)]; subplot(221); plt = @(x) (((x))); pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_e^{00}|$'); subplot(222); plt = @(x) ((x)); pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) ((x)); pclr = pcolor((RR),(ZZ),plt(phi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat\phi|$'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it5D))); legend('$|S_i^{00}|$'); end save_figure end %% t0 = 0; skip_ = 1; DELAY = 0.01*skip_; FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D); %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% Density ion GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$'; create_gif end if 0 %% Density electron GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$'; create_gif end if 0 %% Phi GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1; FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$'; create_gif end if 0 %% Density ion frequency GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$'; create_gif end if 0 %% Density ion frequency @ kr = 0 GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.'; X = (kz); T = Ts2D; YMIN = -.1; YMAX = 1.1; XMIN = min(kz); XMAX = max(kz); FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r$'; create_gif_1D % %% PJ ion moment frequency % p_ = 0+1; j_ = 3+1; % GIFNAME = ['Ni',num2str(p_),num2str(j_),sprintf('_%.2d',JOBNUM)]; INTERP = 0; % FIELD =ifftshift(abs(squeeze(Nipj(p_,j_,:,:,:))),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts; % FIELDNAME = ['$N_i^{',num2str(p_-1),num2str(j_-1),'}$']; XNAME = '$k_r$'; YNAME = '$k_z$'; % create_gif % %% non linear term frequ. space % GIFNAME = ['Si00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; % FIELD = fftshift((abs(Si00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts; % FIELDNAME = '$S_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$'; % create_gif % %% Electron mode growth % GIFNAME = ['norm_Nepj',sprintf('_%.2d',JOBNUM)]; INTERP = 1; % FIELD = Ne_norm; X = PE; Y = JE; T = Ts; % FIELDNAME = '$\max_k{\gamma_e}$'; XNAME = '$p$'; YNAME = '$j$'; % create_gif end %% if 0 %% Asymmetry error on kr = 0 up = 2:Nkz/2; down = Nkz:-1:Nkz/2+2; eps_asym_r = Ts5D; eps_asym_i = Ts5D; for it = 1:numel(Ts5D) eps_asym_r(it) = max(squeeze(abs(real(Ni00(1,up,it))-real(Ni00(1,down,it))))); eps_asym_i(it) = max(abs(imag(Ni00(1,up,it))+imag(Ni00(1,down,it)))); end figure it = 2; subplot(311) plot(kz(up),real(Ni00(1,up,it))); hold on plot(kz(up),real(Ni00(1,down,it)),'--') ylabel('Real'); title(['$N_i^{00}(t=',num2str(Ts5D(it)),')$']) subplot(312) plot(kz(up),imag(Ni00(1,up,it))); hold on plot(kz(up),-imag(Ni00(1,down,it)),'--') ylabel('Imag'); xlabel('$kz,kr=0$') subplot(313) plot(kz(2:Nkz/2),real(Ni00(1,up,it))-real(Ni00(1,down,it))); hold on plot(kz(2:Nkz/2),imag(Ni00(1,up,it))+imag(Ni00(1,down,it)),'--') legend('Re','Im'); ylabel('err'); xlabel('$kz,kr=0$') %% figure semilogy(Ts5D,eps_asym_r); hold on semilogy(Ts5D,eps_asym_i,'--'); title('$\max_{k_z}|N_i^{00}(0,k_z)-N_i^{00}(0,-k_z)|$'); legend('Re','Im'); grid on; xlabel('$t$') end if 0 %% Growth rate fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)]; [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz/kr(ikr0KH),gkr0kz_Ni00/(kr(ikr0KH)*A0KH),'DisplayName',['J = ',num2str(JMAXI)]) xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$'); xlim([0. 1.0]); ylim([0. 0.6]); save_figure end %% if 0 %% Mode time evolution [~,ik00] = min(abs(kz)); [~,idk] = min(abs(kz-dkz)); [~,ik50] = min(abs(kz-0.5*KR0KH)); [~,ik75] = min(abs(kz-0.7*KR0KH)); [~,ik10] = min(abs(kz-1.00*KR0KH)); plt = @(x) abs(squeeze(x)); fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)]; semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik00)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(idk)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik50)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik75)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik10)/KR0KH)]); hold on xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show'); semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts2D-Ts2D(end))),'k--') plot(tstart*[1 1],ylim,'k-','LineWidth',0.5); plot(tend*[1 1],ylim,'k-','LineWidth',0.5); save_figure end %% if 0 %% Show frame in kspace tf = 20; [~,it] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['krkz_frame',sprintf('_%.2d',JOBNUM)]; subplot(221); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat\phi|$'); subplot(222); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_e^{00}|$'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$'); end save_figure end diff --git a/wk/fort.90 b/wk/fort.90 index c849db7..f684ea3 100644 --- a/wk/fort.90 +++ b/wk/fort.90 @@ -1,63 +1,63 @@ &BASIC nrun = 100000000 dt = 0.01 - tmax = 200 + tmax = 300 RESTART = .false. / &GRID pmaxe =2 jmaxe = 1 pmaxi = 2 jmaxi = 1 Nr = 256 Lr = 66 Nz = 256 Lz = 66 kpar = 0 / &OUTPUT_PAR nsave_0d = 0 nsave_1d = 0 nsave_2d = 50 nsave_5d = 1000 write_Na00 = .true. write_moments = .true. write_phi = .true. write_non_lin = .true. write_doubleprecision = .true. - resfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/outputs' - rstfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint' - job2load = 1 + resfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/outputs' + rstfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint' + job2load = 0 / &MODEL_PAR ! Collisionality CO = -2 DK = .false. NON_LIN = .true. mu = 0.0005 nu = 0.01 tau_e = 1 tau_i = 1 sigma_e = 0.023338 sigma_i = 1 q_e = -1 q_i = 1 eta_n = 1 eta_T = 0 - eta_B = 0.5 + eta_B = 0.4 lambdaD = 0 kr0KH = 0 A0KH = 0 / &INITIAL_CON only_Na00 =.false. initback_moments =0 initnoise_moments =1e-05 iseed =42 selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5' eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' / &TIME_INTEGRATION_PAR numerical_scheme='RK4' / \ No newline at end of file diff --git a/wk/load_results.m b/wk/load_results.m index 5fccf08..43bdbde 100644 --- a/wk/load_results.m +++ b/wk/load_results.m @@ -1,21 +1,22 @@ %% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); disp(['Loading ',filename]) % Loading from output file -% CPUTIME = h5readatt(filename,'/data/input','cpu_time'); +CPUTIME = h5readatt(filename,'/data/input','cpu_time'); +DT_SIM = h5readatt(filename,'/data/input','dt'); if strcmp(OUTPUTS.write_moments,'.true.') [Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i'); [Nepj, Pe, Je ] = load_5D_data(filename, 'moments_e'); end [Ni00, kr, kz, Ts2D, dt2D] = load_2D_data(filename, 'Ni00'); Ne00 = load_2D_data(filename, 'Ne00'); %Pi = [0]; Ji = Pi; Pe = Pi; Je = Pi; % Nipj = zeros(1,1,numel(kr),numel(kz),numel(Ts5D)); % Nepj = Nipj; % Nipj(1,1,:,:,:) = Ni00; Nepj(1,1,:,:,:) = Ne00; PHI = load_2D_data(filename, 'phi'); if strcmp(OUTPUTS.write_non_lin,'.true.') Sipj = load_5D_data(filename, 'Sipj'); Sepj = load_5D_data(filename, 'Sepj'); end \ No newline at end of file diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m index 4033fbd..4339223 100644 --- a/wk/parameters_ZP.m +++ b/wk/parameters_ZP.m @@ -1,43 +1,42 @@ clear all; addpath(genpath('../matlab')) % ... add -default_plots_options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 1e-2; % Collision frequency TAU = 1.0; % e/i temperature ratio -ETAB = 0.5; % Magnetic gradient +ETAB = 0.4; % Magnetic gradient ETAN = 1.0; % Density gradient ETAT = 0.0; % Temperature gradient MU = 5e-4; % Hyper diffusivity coefficient NOISE0 = 1.0e-5; %% GRID PARAMETERS N = 256; % Frequency gridpoints (Nkr = N/2) L = 66; % Size of the squared frequency domain PMAXE = 2; % Highest electron Hermite polynomial degree JMAXE = 1; % Highest '' Laguerre '' PMAXI = 2; % Highest ion Hermite polynomial degree JMAXI = 1; % Highest '' Laguerre '' %% TIME PARAMETERS -TMAX = 200; % Maximal time unit +TMAX = 300; % Maximal time unit DT = 1e-2; % Time step SPS2D = 2; % Sampling per time unit for 2D arrays SPS5D = 0.1; % Sampling per time unit for 5D arrays RESTART = 0; % To restart from last checkpoint -JOB2LOAD= 01; +JOB2LOAD= 0; %% OPTIONS SIMID = 'ZP'; % Name of the simulation CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. NO_E = 0; % Remove electrons dynamic % DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1) KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) setup