diff --git a/matlab/COSOlver_matrices_analysis.m b/matlab/COSOlver_matrices_analysis.m index 1514d48..083eb5e 100644 --- a/matlab/COSOlver_matrices_analysis.m +++ b/matlab/COSOlver_matrices_analysis.m @@ -1,72 +1,72 @@ addpath(genpath('../matlab')) % ... add %% Grid configuration -N = 10; % Frequency gridpoints (Nkr = N/2) +N = 10; % Frequency gridpoints (Nkx = N/2) L = 120; % Size of the squared frequency domain dk = 2*pi/L; kmax = N/2*dk; -kr = dk*(0:N/2); -kz = dk*(0:N/2); -[KZ, KR]= meshgrid(kz,kr); -KPERP = sqrt(KR.^2 + KZ.^2); -kperp = reshape(KPERP,[1,numel(kr)^2]); +kx = dk*(0:N/2); +ky = dk*(0:N/2); +[ky, kx]= meshgrid(ky,kx); +KPERP = sqrt(kx.^2 + ky.^2); +kperp = reshape(KPERP,[1,numel(kx)^2]); kperp = uniquetol(kperp,1e-14); Nperp = numel(kperp); COSOlver_path = '../../Documents/MoliSolver/COSOlver/'; COSOLVER.pmaxe = 10; COSOLVER.jmaxe = 5; COSOLVER.pmaxi = 10; COSOLVER.jmaxi = 5; COSOLVER.neFLR = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2)); COSOLVER.neFLRs = 0; % ... only for GK abel COSOLVER.npeFLR = 0; % ... only for GK abel COSOLVER.niFLRs = 0; % ... only for GK abel COSOLVER.npiFLR = 0; % ... only for GK abel COSOLVER.gk = 1; COSOLVER.co = 3; if 0 %% plot the kperp distribution figure plot(kperp) end %% Load the matrices C_self_i = zeros((COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),(COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),Nperp); for n_ = 1:Nperp COSOLVER.kperp = kperp(n_); k_string = sprintf('%0.4f',kperp(n_)); mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),... '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),... '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),... '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),... '_JE_12',... '_NFLR_',num2str(COSOLVER.neFLR),... '_kperp_',k_string,'.h5']; tmp = h5read(mat_file_name,'/Caapj/Ceepj'); C_self_i(:,:,n_) = tmp; end %% Post processing dC_self_i = diff(C_self_i,1,3); gvar_dC = squeeze(sum(sum(abs(dC_self_i),2),1)); %% Plots %% all coeffs evolution figure for ip_ = 1:COSOLVER.pmaxi+1 for ij_ = 1:COSOLVER.jmaxi+1 plot(kperp,squeeze(C_self_i(ip_,ij_,:)),'o'); hold on; end end %% global matrix variation figure; kperpmid = 0.5*(kperp(2:end)+kperp(1:end-1)); plot(kperpmid,gvar_dC./diff(kperp)'); diff --git a/matlab/check_checkpoint.m b/matlab/check_checkpoint.m deleted file mode 100644 index 8265cf4..0000000 --- a/matlab/check_checkpoint.m +++ /dev/null @@ -1,19 +0,0 @@ -function [] = check_checkpoint(cp_ID, OUTPUTS) - -cp_fname = OUTPUTS.rstfile0; -cp_fname = [cp_fname(2:end-1),'_%.2d.h5']; -cp_fname = sprintf(cp_fname,cp_ID) - -tmp = h5read(cp_fname,'/Basic/moments_i'); -Nipj_cp = tmp.real + 1i * tmp.imaginary; -t_cp = h5readatt(cp_fname,'/Basic','time'); - -[~,~,Nr,Nz] = size(Nipj_cp); - -% ni00_cp = fftshift(ifft2(half_2_full_cc_2D(squeeze(Nipj_cp(1,1,:,:))),'symmetric')); -ni00_cp = real(fftshift(ifft2(squeeze(Nipj_cp(1,1,:,:)),Nz,Nz))); -fig = figure; - pclr=pcolor(transpose(ni00_cp(:,:)));set(pclr, 'edgecolor','none'); colorbar; - title(['$n_i^{00}$',', $t \approx$', sprintf('%.3d',ceil(t_cp))]); - xlabel('$i_r$'); ylabel('$i_z$'); -end \ No newline at end of file diff --git a/matlab/compile_results.m b/matlab/compile_results.m index eac0efb..b79637a 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,139 +1,139 @@ CONTINUE = 1; -JOBNUM = 0; JOBFOUND = 0; +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= []; % L_EVOL = []; % DT_EVOL = []; % % FIELDS Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; GGAMMA_ = []; PGAMMA_ = []; PHI_ = []; DENS_E_ = []; DENS_I_ = []; TEMP_E_ = []; TEMP_I_ = []; Ts0D_ = []; -Ts2D_ = []; +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 - Ts2D_ = cat(1,Ts2D_,Ts2D); + Ts3D_ = cat(1,Ts3D_,Ts3D); end if W_PHI - PHI_ = cat(3,PHI_,PHI); + PHI_ = cat(4,PHI_,PHI); end if W_NA00 - Ni00_ = cat(3,Ni00_,Ni00); - Ne00_ = cat(3,Ne00_,Ne00); + Ni00_ = cat(4,Ni00_,Ni00); + Ne00_ = cat(4,Ne00_,Ne00); end if W_DENS - DENS_E_ = cat(3,DENS_E_,DENS_E); - DENS_I_ = cat(3,DENS_I_,DENS_I); + DENS_E_ = cat(4,DENS_E_,DENS_E); + DENS_I_ = cat(4,DENS_I_,DENS_I); end if W_TEMP - TEMP_E_ = cat(3,TEMP_E_,TEMP_E); - TEMP_I_ = cat(3,TEMP_I_,TEMP_I); + 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(5,Nipj_,Nipj); - Nepj_ = cat(5,Nepj_,Nepj); + Nipj_ = cat(6,Nipj_,Nipj); + Nepj_ = cat(6,Nepj_,Nepj); end if W_SAPJ - Sipj_ = cat(5,Sipj_,Sipj); - Sepj_ = cat(5,Sepj_,Sepj); + 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]; 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_; Ts2D = Ts2D_; +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 filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); \ No newline at end of file diff --git a/matlab/compute_Sapj.m b/matlab/compute_Sapj.m deleted file mode 100644 index 0a708aa..0000000 --- a/matlab/compute_Sapj.m +++ /dev/null @@ -1,37 +0,0 @@ -function [ Sapj ] = compute_Sapj(p, j, Kr, Kz, Napj, specie, phi, MODEL, GRID) -%COMPUT_SAPJ compute the non linear term for moment pj specie a -% perform a convolution by product of inverse fourier transform and then -% put the results back into k space with an FFT - Jmax = GRID.jmaxe * (specie=='e')... - + GRID.jmaxi * (specie=='i'); - %padding - Pad = 2.0; - Sapj = zeros(numel(Kr), numel(Kz)); - F = zeros(numel(Kr), numel(Kz)); - G = zeros(numel(Kr), numel(Kz)); - for n = 0:Jmax % Sum over Laguerre - - for ikr = 1:numel(Kr) - for ikz = 1:numel(Kz) - kr = Kr(ikr); kz = Kz(ikz); - BA = sqrt(kr^2+kz^2)*... - (MODEL.sigma_e*sqrt(2) * (specie=='e')... - + sqrt(2*MODEL.tau_i) * (specie=='i')); - - %First conv term - F(ikr,ikz) = (kz-kr).*phi(ikr,ikz).*kernel(n,BA); - %Second conv term - G(ikr,ikz) = 0.0; - for s = 0:min(n+j,Jmax) - G(ikr,ikz) = ... - G(ikr,ikz) + dnjs(n,j,s) .* squeeze(Napj(p+1,s+1,ikr,ikz)); - end - G(ikr,ikz) = (kz-kr) .* G(ikr,ikz); - end - end - %Conv theorem - Sapj = Sapj + conv_thm_2D(F,G,Pad); - end -end - - diff --git a/matlab/conv_thm_2D.m b/matlab/conv_thm_2D.m deleted file mode 100644 index 6f2b445..0000000 --- a/matlab/conv_thm_2D.m +++ /dev/null @@ -1,14 +0,0 @@ -function [ conv ] = conv_thm_2D( F, G, Pad ) -%conv_thm_2D computes the convolution between F and G with a zero pading Pad -% kspace -> real -> product -> kspace - -[Nr, Nz] = size(F); - -f = ifft2(F,Nr*Pad, Nz*Pad); -g = ifft2(G,Nr*Pad, Nz*Pad); - -conv_pad = fft2(f.*g); % convolution becomes product - -conv = conv_pad(1:Nr,1:Nz); % remove padding -end - diff --git a/matlab/create_gif.m b/matlab/create_gif.m index ca3e37b..cd020f8 100644 --- a/matlab/create_gif.m +++ b/matlab/create_gif.m @@ -1,56 +1,60 @@ title1 = GIFNAME; FIGDIR = BASIC.RESDIR; GIFNAME = [FIGDIR, GIFNAME,'.gif']; +XNAME = latexize(XNAME); +YNAME = latexize(YNAME); +FIELDNAME = latexize(FIELDNAME); + % 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(bluewhitered) 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([-1,1]); 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_mov.m b/matlab/create_mov.m index cebdb62..ca28c3c 100644 --- a/matlab/create_mov.m +++ b/matlab/create_mov.m @@ -1,52 +1,56 @@ title1 = GIFNAME; FIGDIR = BASIC.RESDIR; GIFNAME = [FIGDIR, GIFNAME,'.mp4']; +XNAME = latexize(XNAME); +YNAME = latexize(YNAME); +FIELDNAME = latexize(FIELDNAME); + vidfile = VideoWriter(GIFNAME,'Uncompressed AVI'); open(vidfile); % 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 figure('Color','white','Position', [100, 100, 400, 400]); pcolor(X,Y,FIELD(:,:,1)); % to set up % colormap gray colormap(bluewhitered) axis tight manual % this ensures that getframe() returns a consistent size if 1 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 1 shading interp; end set(pclr, 'edgecolor','none'); axis square; % caxis([-max(max(abs(FIELD(:,:,n)))),max(max(abs(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(gcf); writeVideo(vidfile,frame); % 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(['Video saved @ : ',GIFNAME]) close(vidfile); end diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m index 56e500d..f8d5496 100644 --- a/matlab/default_plots_options.m +++ b/matlab/default_plots_options.m @@ -1,10 +1,11 @@ %% 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 = {'-','-.','--',':'}; -marker_styles = {'^','v','o','s'}; \ No newline at end of file +marker_styles = {'^','v','o','s'}; +latexize = @(x) ['$',x,'$']; diff --git a/matlab/fort.90 b/matlab/fort.90 deleted file mode 100644 index 321aad3..0000000 --- a/matlab/fort.90 +++ /dev/null @@ -1,6 +0,0 @@ -&GRID - nr = 64 - Lr = 10 - nz = 64 - Lz = 10 -/ diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m index 47ef91a..b013a63 100644 --- a/matlab/load_2D_data.m +++ b/matlab/load_2D_data.m @@ -1,16 +1,16 @@ function [ data, time, dt ] = load_2D_data( filename, variablename ) %LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ time = h5read(filename,'/data/var2d/time'); - kr = h5read(filename,'/data/grid/coordkr'); - kz = h5read(filename,'/data/grid/coordkz'); + kx = h5read(filename,'/data/grid/coordkx'); + ky = h5read(filename,'/data/grid/coordky'); dt = h5readatt(filename,'/data/input','dt'); cstart= h5readatt(filename,'/data/input','start_iframe2d'); - data = zeros(numel(kr),numel(kz),numel(time)); + data = zeros(numel(kx),numel(ky),numel(time)); for it = 1:numel(time) tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]); data(:,:,it) = tmp.real + 1i * tmp.imaginary; end end diff --git a/matlab/load_3D_data.m b/matlab/load_3D_data.m new file mode 100644 index 0000000..2921908 --- /dev/null +++ b/matlab/load_3D_data.m @@ -0,0 +1,17 @@ +function [ data, time, dt ] = load_3D_data( filename, variablename ) +%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ + time = h5read(filename,'/data/var3d/time'); + kx = h5read(filename,'/data/grid/coordkx'); + ky = h5read(filename,'/data/grid/coordky'); + z = h5read(filename,'/data/grid/coordz'); + dt = h5readatt(filename,'/data/input','dt'); + cstart= h5readatt(filename,'/data/input','start_iframe3d'); + + data = zeros(numel(kx),numel(ky),numel(z),numel(time)); + for it = 1:numel(time) + tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]); + data(:,:,:,it) = tmp.real + 1i * tmp.imaginary; + end + +end + diff --git a/matlab/load_5D_data.m b/matlab/load_4D_data.m similarity index 70% copy from matlab/load_5D_data.m copy to matlab/load_4D_data.m index e36d767..5b7714e 100644 --- a/matlab/load_5D_data.m +++ b/matlab/load_4D_data.m @@ -1,23 +1,23 @@ -function [ data, time, dt ] = load_5D_data( filename, variablename ) +function [ data, time, dt ] = load_4D_data( filename, variablename ) %LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ time = h5read(filename,'/data/var5d/time'); if strcmp(variablename,'moments_e') || strcmp(variablename,'Sepj') p = h5read(filename,'/data/grid/coordp_e'); j = h5read(filename,'/data/grid/coordj_e'); else p = h5read(filename,'/data/grid/coordp_i'); j = h5read(filename,'/data/grid/coordj_i'); end - kr = h5read(filename,'/data/grid/coordkr'); - kz = h5read(filename,'/data/grid/coordkz'); + kx = h5read(filename,'/data/grid/coordkx'); + ky = h5read(filename,'/data/grid/coordky'); dt = h5readatt(filename,'/data/input','dt'); cstart= h5readatt(filename,'/data/input','start_iframe5d'); - data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time)); + data = zeros(numel(p),numel(j),numel(kx),numel(ky),numel(time)); for it = 1:numel(time) tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]); - data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; + data(:,:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; end end \ No newline at end of file diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m index e36d767..2de0553 100644 --- a/matlab/load_5D_data.m +++ b/matlab/load_5D_data.m @@ -1,23 +1,24 @@ function [ data, time, dt ] = load_5D_data( filename, variablename ) %LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ time = h5read(filename,'/data/var5d/time'); if strcmp(variablename,'moments_e') || strcmp(variablename,'Sepj') p = h5read(filename,'/data/grid/coordp_e'); j = h5read(filename,'/data/grid/coordj_e'); else p = h5read(filename,'/data/grid/coordp_i'); j = h5read(filename,'/data/grid/coordj_i'); end - kr = h5read(filename,'/data/grid/coordkr'); - kz = h5read(filename,'/data/grid/coordkz'); + kx = h5read(filename,'/data/grid/coordkx'); + ky = h5read(filename,'/data/grid/coordky'); + z = h5read(filename,'/data/grid/coordz'); dt = h5readatt(filename,'/data/input','dt'); cstart= h5readatt(filename,'/data/input','start_iframe5d'); - data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time)); + data = zeros(numel(p),numel(j),numel(kx),numel(ky),numel(z),numel(time)); for it = 1:numel(time) tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]); - data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; + data(:,:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; end end \ No newline at end of file diff --git a/matlab/load_grid_data.m b/matlab/load_grid_data.m index 4318157..be96424 100644 --- a/matlab/load_grid_data.m +++ b/matlab/load_grid_data.m @@ -1,9 +1,10 @@ -function [ pe, je, pi, ji, kr, kz ] = load_grid_data( filename ) +function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data_3D( filename ) %LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ pe = h5read(filename,'/data/grid/coordp_e'); je = h5read(filename,'/data/grid/coordj_e'); pi = h5read(filename,'/data/grid/coordp_i'); ji = h5read(filename,'/data/grid/coordj_i'); - kr = h5read(filename,'/data/grid/coordkr'); - kz = h5read(filename,'/data/grid/coordkz'); + kx = h5read(filename,'/data/grid/coordkx'); + ky = h5read(filename,'/data/grid/coordky'); + z = h5read(filename,'/data/grid/coordz'); end \ No newline at end of file diff --git a/matlab/load_params.m b/matlab/load_params.m index 28ce1ff..336c305 100644 --- a/matlab/load_params.m +++ b/matlab/load_params.m @@ -1,56 +1,56 @@ CO = h5readatt(filename,'/data/input','CO'); ETAB = h5readatt(filename,'/data/input','eta_B'); ETAN = h5readatt(filename,'/data/input','eta_n'); ETAT = h5readatt(filename,'/data/input','eta_T'); PMAXI = h5readatt(filename,'/data/input','pmaxi'); JMAXI = h5readatt(filename,'/data/input','jmaxi'); PMAXE = h5readatt(filename,'/data/input','pmaxe'); JMAXE = h5readatt(filename,'/data/input','jmaxe'); NON_LIN = h5readatt(filename,'/data/input','NON_LIN'); NU = h5readatt(filename,'/data/input','nu'); -NR = h5readatt(filename,'/data/input','nr'); -NZ = h5readatt(filename,'/data/input','nz'); -L = h5readatt(filename,'/data/input','Lr'); +Nx = h5readatt(filename,'/data/input','Nx'); +Ny = h5readatt(filename,'/data/input','Ny'); +L = h5readatt(filename,'/data/input','Lx'); CLOS = h5readatt(filename,'/data/input','CLOS'); DT_SIM = h5readatt(filename,'/data/input','dt'); MU = h5readatt(filename,'/data/input','mu'); % MU = str2num(filename(end-18:end-14)); %bad... W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y'; W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y'; W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y'; W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y'; W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y'; if NON_LIN == 'y' NON_LIN = 1; else NON_LIN = 0; end if (CO == -3); CONAME = 'PADK'; elseif(CO == -2); CONAME = 'SGDK'; elseif(CO == -1); CONAME = 'DGDK'; elseif(CO == 0); CONAME = 'LB'; elseif(CO == 1); CONAME = 'DGGK'; elseif(CO == 2); CONAME = 'SGGK'; elseif(CO == 3); CONAME = 'PAGK'; 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_',... +degngrad = [degngrad,'_eta_%1.1f_nu_%0.0e_',... CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; -degngrad = sprintf(degngrad,[NU,MU]); +degngrad = sprintf(degngrad,[ETAB/ETAN,NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end -resolution = [num2str(NR),'x',num2str(NZ/2),'_']; +resolution = [num2str(Nx),'x',num2str(Ny/2),'_']; gridname = ['L_',num2str(L),'_']; PARAMS = [resolution,gridname,degngrad]; % BASIC.RESDIR = [SIMDIR,PARAMS,'/']; diff --git a/matlab/load_results.m b/matlab/load_results.m index 412b92a..fd7e85d 100644 --- a/matlab/load_results.m +++ b/matlab/load_results.m @@ -1,51 +1,51 @@ %% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp(['Loading ',filename]) % Loading from output file CPUTIME = h5readatt(filename,'/data/input','cpu_time'); DT_SIM = h5readatt(filename,'/data/input','dt'); -[Pe, Je, Pi, Ji, kr, kz] = load_grid_data(filename); +[Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename); W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y'); W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi') ,'y'); W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00') ,'y'); W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj') ,'y'); W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj') ,'y'); W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens') ,'y'); W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp') ,'y'); if W_GAMMA [ GGAMMA_RI, Ts0D, dt0D] = load_0D_data(filename, 'gflux_ri'); PGAMMA_RI = load_0D_data(filename, 'pflux_ri'); end if W_PHI - [ PHI, Ts2D, dt2D] = load_2D_data(filename, 'phi'); + [ PHI, Ts3D, dt3D] = load_3D_data(filename, 'phi'); end if W_NA00 - [Ni00, Ts2D, dt2D] = load_2D_data(filename, 'Ni00'); - Ne00 = load_2D_data(filename, 'Ne00'); + [Ni00, Ts3D, dt3D] = load_3D_data(filename, 'Ni00'); + Ne00 = load_3D_data(filename, 'Ne00'); end if W_NAPJ [Nipj, Ts5D, dt5D] = load_5D_data(filename, 'moments_i'); [Nepj ] = load_5D_data(filename, 'moments_e'); end if W_SAPJ [Sipj, Ts5D, dt5D] = load_5D_data(filename, 'Sipj'); Sepj = load_5D_data(filename, 'Sepj'); end if W_DENS - [DENS_E, Ts2D, dt2D] = load_2D_data(filename, 'dens_e'); - [DENS_I, Ts2D, dt2D] = load_2D_data(filename, 'dens_i'); + [DENS_E, Ts3D, dt3D] = load_3D_data(filename, 'dens_e'); + [DENS_I, Ts3D, dt3D] = load_3D_data(filename, 'dens_i'); end if W_TEMP - [TEMP_E, Ts2D, dt2D] = load_2D_data(filename, 'temp_e'); - [TEMP_I, Ts2D, dt2D] = load_2D_data(filename, 'temp_i'); + [TEMP_E, Ts3D, dt3D] = load_3D_data(filename, 'temp_e'); + [TEMP_I, Ts3D, dt3D] = load_3D_data(filename, 'temp_i'); end \ No newline at end of file diff --git a/matlab/plot_kernels.m b/matlab/plot_kernels.m index f772984..f7eede8 100644 --- a/matlab/plot_kernels.m +++ b/matlab/plot_kernels.m @@ -1,46 +1,46 @@ %% Kernels kmax=7; nmax=6; -kr_ = linspace(0,kmax,100); +kx_ = linspace(0,kmax,100); figure for n_ = 0:nmax - plot(kr_,kernel(n_,kr_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on; + plot(kx_,kernel(n_,kx_),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on; end ylim_ = ylim; -plot(kr_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); -plot(kr_,J0,'-r','DisplayName','$J_0$'); +plot(kx_(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); +plot(kx_,J0,'-r','DisplayName','$J_0$'); legend('show') %% Bessels and approx vperp = linspace(0,1.5,4); nmax1=5; nmax2=10; kmax=7; figure for i = 1:4 subplot(2,2,i) v_ = vperp(i); - kr_ = linspace(0,kmax,100); + kx_ = linspace(0,kmax,100); - J0 = besselj(0,kr_*v_); - A1 = 1 - kr_.^2*v_^2/4; - K1 = zeros(size(kr_)); - K2 = zeros(size(kr_)); + 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_,kr_).*polyval(LaguerrePoly(n_),v_^2); + K1 = K1 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2); end for n_ = 0:nmax2 - K2 = K2 + kernel(n_,kr_).*polyval(LaguerrePoly(n_),v_^2); + K2 = K2 + kernel(n_,kx_).*polyval(LaguerrePoly(n_),v_^2); end - plot(kr_,J0,'-k','DisplayName','$J_0(kv)$'); hold on; - plot(kr_,A1,'-r','DisplayName','$1 - k^2 v^2/4$'); - plot(kr_,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']); - plot(kr_,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']); + 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(kr_(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 diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plots/plot_kperp_spectrum.m new file mode 100644 index 0000000..2370dd0 --- /dev/null +++ b/matlab/plots/plot_kperp_spectrum.m @@ -0,0 +1,37 @@ +[~,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]); +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); +%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]) +% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on; +plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on; +plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on; +plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on; +plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$'); +plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$'); +set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on +xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest') +title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... +', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... +' $\mu_{hd}=$',num2str(MU)]}); +xlim([0.1,10]); +save_figure +clear Z_rth phi_kp ni_kp Ti_kp \ No newline at end of file diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m new file mode 100644 index 0000000..ba86525 --- /dev/null +++ b/matlab/plots/plot_radial_transport_and_shear.m @@ -0,0 +1,44 @@ +%Compute steady radial transport +tend = Ts0D(end); tstart = tend - TAVG; +[~,its0D] = min(abs(Ts0D-tstart)); +[~,ite0D] = min(abs(Ts0D-tend)); +SCALE = (2*pi/Nx/Ny)^2; +gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE; +gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE; +% Compute steady shearing rate +tend = Ts3D(end); tstart = tend - TAVG; +[~,its2D] = min(abs(Ts3D-tstart)); +[~,ite2D] = min(abs(Ts3D-tend)); +shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1)); +shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1)); +% plots +fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) + subplot(311) +% yyaxis left + plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on; + plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',... + 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]); + grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$') + ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]); + 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)]); + % + subplot(312) + clr = line_colors(1,:); + lstyle = line_styles(1); +% plt = @(x_) mean(x_,1); + plt = @(x_) x_(1,:); + plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on; + plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on; + plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on; + plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on; + plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',... + 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]); + ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]); + grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show'); + subplot(313) + [TY,TX] = meshgrid(x,Ts3D); + pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar; + caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256)) +save_figure \ No newline at end of file diff --git a/matlab/plots/plot_shear_and_reynold_stress.m b/matlab/plots/plot_shear_and_reynold_stress.m new file mode 100644 index 0000000..e951208 --- /dev/null +++ b/matlab/plots/plot_shear_and_reynold_stress.m @@ -0,0 +1,10 @@ +% trange = 100:200; +figure; +plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3)))); +plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on; +plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on; +% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on; +plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on; +plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on; +% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on; +xlim([-L/2,L/2]); xlabel('$x/\rho_s$'); grid on; legend('show') \ No newline at end of file diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m new file mode 100644 index 0000000..3e263a0 --- /dev/null +++ b/matlab/plots/plot_space_time_diagrams.m @@ -0,0 +1,24 @@ +[~,itstart] = min(abs(Ts3D-tstart)); +[~,itend] = min(abs(Ts3D-tend)); +trange = itstart:itend; +[TY,TX] = meshgrid(x,Ts3D(trange)); +fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) + subplot(211) +% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; + pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; + shading interp + colormap hot; + caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]); + caxis([0.0,cmax]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z'; + xticks([]); ylabel('$x/\rho_s$') +% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$') + 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)]); + subplot(212) + pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; + fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2))); + caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z'; + xlabel('$t c_s/R$'), ylabel('$x/\rho_s$') +% legend('Zonal flow $\langle \partial_r\phi\rangle_z$') +save_figure \ No newline at end of file diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plots/plot_time_evolution_and_gr.m new file mode 100644 index 0000000..e016520 --- /dev/null +++ b/matlab/plots/plot_time_evolution_and_gr.m @@ -0,0 +1,64 @@ +fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS]; +set(gcf, 'Position', [100, 100, 900, 800]) +subplot(111); + suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... + ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... + ' $\mu_{hd}=$',num2str(MU)]); + subplot(421); + for ip = 1:Pe_max + for ij = 1:Je_max + plt = @(x) squeeze(x(ip,ij,:)); + plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$']; + clr = line_colors(min(ip,numel(line_colors(:,1))),:); + lstyle = line_styles(min(ij,numel(line_styles))); + semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,... + 'Color',clr,'LineStyle',lstyle{1}); hold on; + end + end + grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); + subplot(423) + for ip = 1:Pi_max + for ij = 1:Ji_max + plt = @(x) squeeze(x(ip,ij,:)); + plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$']; + clr = line_colors(min(ip,numel(line_colors(:,1))),:); + lstyle = line_styles(min(ij,numel(line_styles))); + semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,... + 'Color',clr,'LineStyle',lstyle{1}); hold on; + end + end + grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$') + subplot(222) + plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on; + plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2); + legend(['Gyro. flux';'Part. flux']); + grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$') + if(~isnan(max(max(g_I(1,:,:))))) + subplot(223) + plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on; + plot(kx,max(g_II(:,1,:),[],3),'x-','DisplayName','Second. instability'); hold on; + plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA'); + grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show'); + ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]); + shearplot = 426; phiplot = 428; + else + shearplot = 223; phiplot = 224; + end + subplot(shearplot) + plt = @(x) mean(x,1); + clr = line_colors(min(ip,numel(line_colors(:,1))),:); + lstyle = line_styles(min(ij,numel(line_styles))); + plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s)$'); hold on; + plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on; + plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on; + plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on; + grid on; xlabel('$t c_s/R$'); ylabel('$shear$'); + subplot(phiplot) + clr = line_colors(min(ip,numel(line_colors(:,1))),:); + lstyle = line_styles(min(ij,numel(line_styles))); + plot(Ts3D,plt(phi_maxx_maxy),'DisplayName','$\max_{r,z}(\phi)$'); hold on; + plot(Ts3D,plt(phi_maxx_avg),'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on; + plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on; + plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on; + grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$'); +save_figure \ No newline at end of file diff --git a/matlab/post_processing.m b/matlab/post_processing.m new file mode 100644 index 0000000..ea25ef8 --- /dev/null +++ b/matlab/post_processing.m @@ -0,0 +1,191 @@ + +%% 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); +Ns3D = numel(Ts3D); +% renaming and reshaping quantity of interest +Ts5D = Ts5D'; +Ts3D = Ts3D'; + +%% Build grids +Nkx = numel(kx); Nky = numel(ky); +[KY,KX] = meshgrid(ky,kx); +Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky); +dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1); +KPERP2 = KY.^2+KX.^2; +[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky)); + +Lk = max(Lkx,Lky); +Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z); +dx = 2*pi/Lk; dy = 2*pi/Lk; dz = q0*2*pi/Nz; +x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x); +y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y); +z = dz * (1:Nz); +[Y_XY,X_XY] = meshgrid(y,x); +[Z_XZ,X_XZ] = meshgrid(z,x); +[Z_YZ,Y_YZ] = meshgrid(z,y); + +%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +disp('Analysis :') +disp('- iFFT') +% IFFT (Lower case = real space, upper case = frequency space) +ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density +ni00 = zeros(Nx,Ny,Nz,Ns3D); +dzTe = zeros(Nx,Ny,Nz,Ns3D); +dzTi = zeros(Nx,Ny,Nz,Ns3D); +dzni = zeros(Nx,Ny,Nz,Ns3D); +np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density +si00 = zeros(Nx,Ny,Nz,Ns5D); +phi = zeros(Nx,Ny,Nz,Ns3D); +dens_e = zeros(Nx,Ny,Nz,Ns3D); +dens_i = zeros(Nx,Ny,Nz,Ns3D); +temp_e = zeros(Nx,Ny,Nz,Ns3D); +temp_i = zeros(Nx,Ny,Nz,Ns3D); +drphi = zeros(Nx,Ny,Nz,Ns3D); +dzphi = zeros(Nx,Ny,Nz,Ns3D); +dr2phi = zeros(Nx,Ny,Nz,Ns3D); + +for it = 1:numel(Ts3D) + for iz = 1:numel(z) + NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); + ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny))); + ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny))); + phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny))); + drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny))); + dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny))); + dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny))); + if(W_DENS && W_TEMP) + DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it); + TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it); + dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny))); + dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny))); + dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny))); + dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny))); + dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny))); + temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny))); + temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny))); + end + end +end + +Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space + +for it = 1:numel(Ts5D) + [~, it2D] = min(abs(Ts3D-Ts5D(it))); + Np_i(:,:,it) = 0; + for ij = 1:Nji + Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1)); + Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it)); + end + + np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny))); +end + +% Post processing +disp('- post processing') +% gyrocenter and particle flux from real space +GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = +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); + +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_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)))); + + 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))))); + + 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 92352d9..517c921 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,135 +1,144 @@ %% ________________________________________________________________________ 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.Nr = N; % r grid resolution -GRID.Lr = L; % r length -GRID.Nz = N * (1-KREQ0) + KREQ0; % z '' -GRID.Lz = L * (1-KREQ0); % z '' -GRID.kpar = KPAR; +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 % 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*KR0KH),'_A_',num2str(A0KH)]; end; +% 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 disp('Warning, FCGK not implemented yet') elseif (CO == -1) % DGDK disp('Warning, DGDK not debugged') end % Naming and creating input file if (CO == -3); CONAME = 'PADK'; elseif(CO == -2); CONAME = 'SGDK'; elseif(CO == -1); CONAME = 'DGDK'; elseif(CO == 0); CONAME = 'LBDK'; elseif(CO == 1); CONAME = 'DGGK'; elseif(CO == 2); CONAME = 'SGGK'; elseif(CO == 3); CONAME = 'PAGK'; 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,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; + CONAME,'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end -resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_']; -gridname = ['L_',num2str(L),'_']; +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),'_q_',num2str(q),'_']; +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/matlab/write_fort90.m b/matlab/write_fort90.m index 19ccdf5..366ee16 100644 --- a/matlab/write_fort90.m +++ b/matlab/write_fort90.m @@ -1,82 +1,86 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC) % Write the input script "fort.90" with desired parameters INPUT = 'fort.90'; fid = fopen(INPUT,'wt'); fprintf(fid,'&BASIC\n'); fprintf(fid,[' nrun = ', num2str(BASIC.nrun),'\n']); fprintf(fid,[' dt = ', num2str(BASIC.dt),'\n']); fprintf(fid,[' tmax = ', num2str(BASIC.tmax),'\n']); fprintf(fid,[' RESTART = ', BASIC.RESTART,'\n']); fprintf(fid,[' maxruntime = ', num2str(BASIC.maxruntime),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&GRID\n'); fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),'\n']); fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']); fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']); fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']); -fprintf(fid,[' Nr = ', num2str(GRID.Nr),'\n']); -fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']); +fprintf(fid,[' Nx = ', num2str(GRID.Nx),'\n']); +fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']); +fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']); +fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']); fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']); -fprintf(fid,[' Lz = ', num2str(GRID.Lz),'\n']); -fprintf(fid,[' kpar = ', num2str(GRID.kpar),'\n']); +fprintf(fid,[' q0 = ', num2str(GRID.q0),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&OUTPUT_PAR\n'); fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']); fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']); fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']); +fprintf(fid,[' nsave_3d = ', num2str(OUTPUTS.nsave_3d),'\n']); fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']); fprintf(fid,[' nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']); fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']); fprintf(fid,[' write_gamma = ', OUTPUTS.write_gamma,'\n']); fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']); fprintf(fid,[' write_Na00 = ', OUTPUTS.write_Na00,'\n']); fprintf(fid,[' write_Napj = ', OUTPUTS.write_Napj,'\n']); fprintf(fid,[' write_Sapj = ', OUTPUTS.write_Sapj,'\n']); fprintf(fid,[' write_dens = ', OUTPUTS.write_dens,'\n']); fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']); fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']); fprintf(fid,[' rstfile0 = ', OUTPUTS.rstfile0,'\n']); fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&MODEL_PAR\n'); fprintf(fid,' ! Collisionality\n'); fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']); fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']); fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']); fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']); fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']); fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']); fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']); fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']); fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']); fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']); fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']); fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']); fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']); fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']); fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),'\n']); fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),'\n']); fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),'\n']); fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&INITIAL_CON\n'); fprintf(fid,[' INIT_NOISY_PHI =', INITIAL.init_noisy_phi,'\n']); fprintf(fid,[' INIT_ZF =', num2str(INITIAL.INIT_ZF),'\n']); +fprintf(fid,[' WIPE_TURB =', INITIAL.wipe_turb,'\n']); +fprintf(fid,[' INIT_BLOB =', INITIAL.init_blob,'\n']); fprintf(fid,[' init_background =', num2str(INITIAL.init_background),'\n']); fprintf(fid,[' init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']); fprintf(fid,[' iseed =', num2str(INITIAL.iseed),'\n']); fprintf(fid,[' mat_file =', INITIAL.mat_file,'\n']); fprintf(fid,'/\n'); fprintf(fid,'&TIME_INTEGRATION_PAR\n'); fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']); fprintf(fid,'/'); fclose(fid); system(['cp fort.90 ',BASIC.RESDIR,'/.']); end diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m index 98c1e35..62a48ac 100644 --- a/matlab/write_sbash_daint.m +++ b/matlab/write_sbash_daint.m @@ -1,57 +1,57 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $SCRATCH/HeLaZ/wk\n',... ... 'cd $SCRATCH/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... 'cp $HOME/HeLaZ/wk/fort.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... 'jid=$(sbatch batch_script.sh)\n',... 'echo $jid\n',... 'echo to check output log :\n',... 'echo tail -f $SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --job-name="',CLUSTER.JNAME,'"\n',... '#SBATCH --time=', CLUSTER.TIME,'\n',... '#SBATCH --nodes=', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',... '#SBATCH --ntasks-per-core=', CLUSTER.NTPC,'\n',... '#SBATCH --mem=', CLUSTER.MEM,'\n',... '#SBATCH --error=err.txt\n',... '#SBATCH --output=out.txt\n',... '#SBATCH --account="s882"\n',... '#SBATCH --constraint=mc\n',... '#SBATCH --hint=nomultithread\n',... '#SBATCH --partition=',CLUSTER.PART,'\n',... 'export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n',... ...% '#SBATCH --job-name=',PARAMS,'\n\n',... 'module purge\n',... 'module load PrgEnv-intel\n',... 'module load cray-hdf5-parallel\n',... 'module load cray-mpich\n',... 'module load craype-x86-skylake\n',... 'module load cray-fftw\n',... -'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]); +'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]); %'srun ./../../../bin/helaz']); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk'); \ No newline at end of file diff --git a/matlab/write_sbash_izar.m b/matlab/write_sbash_izar.m index 01c7334..02224f4 100644 --- a/matlab/write_sbash_izar.m +++ b/matlab/write_sbash_izar.m @@ -1,46 +1,46 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $SCRATCH/ahoffman/HeLaZ/wk\n',... ... 'cd $SCRATCH/ahoffman/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... 'cp $HOME/HeLaZ/wk/fort.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... 'sbatch batch_script.sh\n',... 'echo tail -f $SCRATCH/ahoffman/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --chdir $SCRATCH/ahoffman/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'\n',... '#SBATCH --job-name ',CLUSTER.JNAME,'\n',... '#SBATCH --time ', CLUSTER.TIME,'\n',... '#SBATCH --nodes ', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task ', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node ', CLUSTER.NTPN,'\n',... '#SBATCH --mem ', CLUSTER.MEM,'\n',... '#SBATCH --partition ',CLUSTER.PART,'\n',... 'module purge\n',... 'module load intel/19.0.5\n',... 'module load intel-mpi/2019.5.281\n',... 'module load fftw/3.3.8-mpi-openmp\n',... 'module load hdf5/1.10.6-mpi\n',... -'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]); +'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@izar.epfl.ch:/home/ahoffman/HeLaZ/wk'); \ No newline at end of file diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index 6fc1d66..397c03b 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,510 +1,145 @@ addpath(genpath('../matlab')) % ... add -for i_ = 1 -% for ETA_ =[0.6:0.1:0.9] -%% Load results +addpath(genpath('../matlab/plots')) % ... add +%% Directory of the simulation if 1% Local results -outfile =''; -outfile ='Blob_diffusion/100x50_L_60_P_2_J_1_eta_Inf_nu_1e-01_DGGK_mu_0e+00'; -% outfile ='test_3D/100x50x10_L_60_q0_1_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_2e-03'; -% outfile ='test_3D/100x50_L_60_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_0e+00'; +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); end if 0% Marconi results outfile =''; outfile =''; outfile =''; outfile =''; outfile =''; outfile =''; 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_3D %Compile the results from first output found to JOBNUMMAX if existing - -%% 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); -Ns3D = numel(Ts3D); -% renaming and reshaping quantity of interest -Ts5D = Ts5D'; -Ts3D = Ts3D'; - -%% Build grids -Nkx = numel(kx); Nky = numel(ky); -[KY,KX] = meshgrid(ky,kx); -Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky); -dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1); -KPERP2 = KY.^2+KX.^2; -[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky)); - -Lk = max(Lkx,Lky); -Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z); -dx = 2*pi/Lk; dy = 2*pi/Lk; dz = q0*2*pi/Nz; -x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x); -y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y); -z = dz * (1:Nz); -[Y_XY,X_XY] = meshgrid(y,x); -[Z_XZ,X_XZ] = meshgrid(z,x); -[Z_YZ,Y_YZ] = meshgrid(z,y); - -%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -disp('Analysis :') -disp('- iFFT') -% IFFT (Lower case = real space, upper case = frequency space) -ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density -ni00 = zeros(Nx,Ny,Nz,Ns3D); -dzTe = zeros(Nx,Ny,Nz,Ns3D); -dzTi = zeros(Nx,Ny,Nz,Ns3D); -dzni = zeros(Nx,Ny,Nz,Ns3D); -np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density -si00 = zeros(Nx,Ny,Nz,Ns5D); -phi = zeros(Nx,Ny,Nz,Ns3D); -dens_e = zeros(Nx,Ny,Nz,Ns3D); -dens_i = zeros(Nx,Ny,Nz,Ns3D); -temp_e = zeros(Nx,Ny,Nz,Ns3D); -temp_i = zeros(Nx,Ny,Nz,Ns3D); -drphi = zeros(Nx,Ny,Nz,Ns3D); -dzphi = zeros(Nx,Ny,Nz,Ns3D); -dr2phi = zeros(Nx,Ny,Nz,Ns3D); - -for it = 1:numel(Ts3D) - for iz = 1:numel(z) - NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); - ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny))); - ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny))); - phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny))); - drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny))); - dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny))); - dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny))); - if(W_DENS && W_TEMP) - DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it); - TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it); - dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny))); - dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny))); - dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny))); - dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny))); - dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny))); - temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny))); - temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny))); - end - end -end - -Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space +compile_results %Compile the results from first output found to JOBNUMMAX if existing -for it = 1:numel(Ts5D) - [~, it2D] = min(abs(Ts3D-Ts5D(it))); - Np_i(:,:,it) = 0; - for ij = 1:Nji - Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1)); - Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it)); - end - - np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny))); -end - -% Post processing -disp('- post processing') -% gyrocenter and particle flux from real space -GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = -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); - -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_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)))); - - 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))))); - - 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)); +%% Post-processing +post_processing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Time evolutions and growth rate -fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS]; -set(gcf, 'Position', [100, 100, 900, 800]) -subplot(111); - suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... - ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... - ' $\mu_{hd}=$',num2str(MU)]); - subplot(421); - for ip = 1:Pe_max - for ij = 1:Je_max - plt = @(x) squeeze(x(ip,ij,:)); - plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$']; - clr = line_colors(min(ip,numel(line_colors(:,1))),:); - lstyle = line_styles(min(ij,numel(line_styles))); - semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,... - 'Color',clr,'LineStyle',lstyle{1}); hold on; - end - end - grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); - subplot(423) - for ip = 1:Pi_max - for ij = 1:Ji_max - plt = @(x) squeeze(x(ip,ij,:)); - plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$']; - clr = line_colors(min(ip,numel(line_colors(:,1))),:); - lstyle = line_styles(min(ij,numel(line_styles))); - semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,... - 'Color',clr,'LineStyle',lstyle{1}); hold on; - end - end - grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$') - subplot(222) - plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on; - plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2); - legend(['Gyro. flux';'Part. flux']); - grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$') - if(~isnan(max(max(g_I(1,:,:))))) - subplot(223) - plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on; - plot(kx,max(g_II(:,1,:),[],3),'x-','DisplayName','Second. instability'); hold on; - plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA'); - grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show'); - ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]); - shearplot = 426; phiplot = 428; - else - shearplot = 223; phiplot = 224; - end - subplot(shearplot) - plt = @(x) mean(x,1); - clr = line_colors(min(ip,numel(line_colors(:,1))),:); - lstyle = line_styles(min(ij,numel(line_styles))); - plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s)$'); hold on; - plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on; - plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on; - plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on; - grid on; xlabel('$t c_s/R$'); ylabel('$shear$'); - subplot(phiplot) - clr = line_colors(min(ip,numel(line_colors(:,1))),:); - lstyle = line_styles(min(ij,numel(line_styles))); - plot(Ts3D,plt(phi_maxx_maxy),'DisplayName','$\max_{r,z}(\phi)$'); hold on; - plot(Ts3D,plt(phi_maxx_avg),'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on; - plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on; - plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on; - grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$'); -save_figure +plot_time_evolution_and_gr end if 1 %% Space time diagramm (fig 11 Ivanov 2020) TAVG = 5000; % Averaging time duration -%Compute steady radial transport -tend = Ts0D(end); tstart = tend - TAVG; -[~,its0D] = min(abs(Ts0D-tstart)); -[~,ite0D] = min(abs(Ts0D-tend)); -SCALE = (2*pi/Nx/Ny)^2; -gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE; -gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE; -% Compute steady shearing rate -tend = Ts3D(end); tstart = tend - TAVG; -[~,its2D] = min(abs(Ts3D-tstart)); -[~,ite2D] = min(abs(Ts3D-tend)); -shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1)); -shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1)); -% plots -fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) - subplot(311) -% yyaxis left - plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on; - plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',... - 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]); - grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$') - ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]); - 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)]); - % - subplot(312) - clr = line_colors(1,:); - lstyle = line_styles(1); -% plt = @(x_) mean(x_,1); - plt = @(x_) x_(1,:); - plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on; - plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on; - plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on; - plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on; - plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',... - 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]); - ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]); - grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show'); - subplot(313) - [TY,TX] = meshgrid(x,Ts3D); - pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar; - caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256)) -save_figure +plot_radial_transport_and_shear end if 0 %% Space time diagramms -tstart = 0; tend = Ts3D(end); -[~,itstart] = min(abs(Ts3D-tstart)); -[~,itend] = min(abs(Ts3D-tend)); -trange = itstart:itend; -[TY,TX] = meshgrid(x,Ts3D(trange)); -fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) - subplot(211) -% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; - pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; - shading interp - colormap hot; - caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]); - caxis([0.0,0.01]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z'; - xticks([]); ylabel('$x/\rho_s$') -% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$') - 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)]); - subplot(212) - pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; - fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2))); - caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z'; - xlabel('$t c_s/R$'), ylabel('$x/\rho_s$') -% legend('Zonal flow $\langle \partial_r\phi\rangle_z$') -save_figure +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; -% trange = 100:200; -figure; -plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3)))); -plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on; -plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on; -% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on; -plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on; -plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on; -% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on; -xlim([-L/2,L/2]); xlabel('$x/\rho_s$'); grid on; legend('show') +plot_shear_and_reynold_stress end if 0 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3) -tstart = 0.8*Ts3D(end); tend = Ts3D(end); -[~,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]); -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); -%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]) -% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on; -plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on; -plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on; -plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on; -plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$'); -plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$'); -set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on -xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest') -title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... -', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... -' $\mu_{hd}=$',num2str(MU)]}); -xlim([0.1,10]); -save_figure -clear Z_rth phi_kp ni_kp Ti_kp +tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window +plot_kperp_spectrum end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options t0 =00; 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 = 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 % 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'; % 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'; 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'; % Chose when to plot it tf = [0 1 2 3]; % Slice 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))),')'] % 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'] % 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 %% -% ZF_fourier_analysis -end \ No newline at end of file +% ZF_fourier_analysis \ No newline at end of file