diff --git a/matlab/compile_results.m b/matlab/compile_results.m index b98ac6c..8beff24 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,268 +1,270 @@ function [DATA] = compile_results(DIRECTORY,JOBNUMMIN,JOBNUMMAX) DATA = {}; CONTINUE = 1; JOBNUM = JOBNUMMIN; JOBFOUND = 0; DATA.TJOB_SE = []; % Start and end times of jobs DATA.NU_EVOL = []; % evolution of parameter nu between jobs DATA.CO_EVOL = []; % evolution of CO DATA.MUx_EVOL = []; % evolution of parameter mu between jobs DATA.MUy_EVOL = []; % evolution of parameter mu between jobs DATA.K_N_EVOL = []; % DATA.L_EVOL = []; % DATA.DT_EVOL = []; % % FIELDS Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; HFLUX_ = []; GGAMMA_ = []; PGAMMA_ = []; PHI_ = []; DENS_E_ = []; DENS_I_ = []; UPAR_E_ = []; UPAR_I_ = []; UPER_E_ = []; UPER_I_ = []; TPAR_E_ = []; TPAR_I_ = []; TPER_E_ = []; TPER_I_ = []; TEMP_E_ = []; TEMP_I_ = []; TEMP_E_ = []; TEMP_I_ = []; Ts0D_ = []; Ts3D_ = []; Ts5D_ = []; Sipj_ = []; Sepj_ = []; Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old; Pi_max=0; Pe_max=0; Ji_max=0; Je_max=0; +DATA.outfilenames = {}; +ii = 1; while(CONTINUE) filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM); if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX) - % Load results of simulation #JOBNUM -% load_results + DATA.outfilenames{ii} = filename; %% 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, kx, ky, z] = load_grid_data(filename); W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y'); W_HF = strcmp(h5readatt(filename,'/data/input','write_hf' ),'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'); KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y'); % 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_RI, Ts0D, ~] = load_0D_data(filename, 'gflux_ri'); PGAMMA_RI = load_0D_data(filename, 'pflux_ri'); GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); clear GGAMMA_RI PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); clear PGAMMA_RI end if W_HF [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x'); HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X end if W_PHI [ PHI, Ts3D, ~] = load_3D_data(filename, 'phi'); PHI_ = cat(4,PHI_,PHI); clear PHI end if W_NA00 if KIN_E Ne00 = load_3D_data(filename, 'Ne00'); Ne00_ = cat(4,Ne00_,Ne00); clear Ne00 end [Ni00, Ts3D, ~] = load_3D_data(filename, 'Ni00'); Ni00_ = cat(4,Ni00_,Ni00); clear Ni00 end if W_DENS if KIN_E [DENS_E, Ts3D, ~] = load_3D_data(filename, 'dens_e'); DENS_E_ = cat(4,DENS_E_,DENS_E); clear DENS_E end [DENS_I, Ts3D, ~] = load_3D_data(filename, 'dens_i'); DENS_I_ = cat(4,DENS_I_,DENS_I); clear DENS_I end if 0 if KIN_E [UPAR_E, Ts3D, ~] = load_3D_data(filename, 'upar_e'); UPAR_E_ = cat(4,UPAR_E_,UPAR_E); clear UPAR_E [UPER_E, Ts3D, ~] = load_3D_data(filename, 'uper_e'); % UPER_E_ = cat(4,UPER_E_,UPER_E); clear UPER_E end [UPAR_I, Ts3D, ~] = load_3D_data(filename, 'upar_i'); UPAR_I_ = cat(4,UPAR_I_,UPAR_I); clear UPAR_I [UPER_I, Ts3D, ~] = load_3D_data(filename, 'uper_i'); UPER_I_ = cat(4,UPER_I_,UPER_I); clear UPER_I end if W_TEMP if KIN_E % [TPAR_E, Ts3D, ~] = load_3D_data(filename, 'Tpar_e'); % TPAR_E_ = cat(4,TPAR_E_,TPAR_E); clear TPAR_E % [TPER_E, Ts3D, ~] = load_3D_data(filename, 'Tper_e'); % TPER_E_ = cat(4,TPER_E_,TPER_E); clear TPER_E [TEMP_E, Ts3D, ~] = load_3D_data(filename, 'temp_e'); TEMP_E_ = cat(4,TEMP_E_,TEMP_E); clear TEMP_E end % [TPAR_I, Ts3D, ~] = load_3D_data(filename, 'Tpar_i'); % TPAR_I_ = cat(4,TPAR_I_,TPAR_I); clear TPAR_I % [TPER_I, Ts3D, ~] = load_3D_data(filename, 'Tper_i'); % TPER_I_ = cat(4,TPER_I_,TPER_I); clear TPER_I [TEMP_I, Ts3D, ~] = load_3D_data(filename, 'temp_i'); TEMP_I_ = cat(4,TEMP_I_,TEMP_I); clear TEMP_I end Ts5D = []; if W_NAPJ [Nipj, Ts5D, ~] = load_5D_data(filename, 'moments_i'); Nipj_ = cat(6,Nipj_,Nipj); clear Nipj if KIN_E Nepj = load_5D_data(filename, 'moments_e'); Nepj_ = cat(6,Nepj_,Nepj); clear Nepj end end if W_SAPJ Sipj_ = cat(6,Sipj_,Sipj); if KIN_E Sepj_ = cat(6,Sepj_,Sepj); end end Ts0D_ = cat(1,Ts0D_,Ts0D); Ts3D_ = cat(1,Ts3D_,Ts3D); Ts5D_ = cat(1,Ts5D_,Ts5D); % Evolution of simulation parameters load_params DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)]; DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU]; DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO]; DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx]; DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy]; DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N]; DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L]; DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM]; - + + ii = ii + 1; 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 %% Build grids Nkx = numel(kx); if Nkx > 1 dkx = kx(2); Lx = 2*pi/dkx; else dkx = 0; Lx = 0; end [~,ikx0] = min(abs(kx)); Nx = 2*Nkx-1; x = linspace(-Lx/2,Lx/2,Nx+1); x = x(1:end-1); Nky = numel(ky); if Nky > 1 dky = ky(2); Ly = 2*pi/dky; else dky = 0; Ly = 0; end [~,iky0] = min(abs(ky)); Ny = Nky; y = linspace(-Ly/2,Ly/2,Ny+1); y = y(1:end-1); Nz = numel(z); [KY,KX] = meshgrid(ky,kx); KPERP2 = KY.^2+KX.^2; %% Add everything in output structure % scaling DATA.scale = (1/Nx/Ny)^2; % Fields DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_; DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_; if(KIN_E) DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_; end DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_; DATA.PHI = PHI_; % grids DATA.Pe = Pe; DATA.Pi = Pi; DATA.Je = Je; DATA.Ji = Ji; DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.x = x; DATA.y = y; DATA.ikx0 = ikx0; DATA.iky0 = iky0; DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky; DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji); DATA.dir = DIRECTORY; DATA.localdir = ['..',DIRECTORY(20:end)]; DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ... ', $\kappa_N=$',num2str(DATA.K_N),', $L=',num2str(DATA.L),'$, $N=',... num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',... num2str(DATA.JMAXI),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),... ',',num2str(DATA.MUy),')']; JOBNUM = LASTJOB; filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM); end \ No newline at end of file diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m index 60fb943..b92ecbe 100644 --- a/matlab/compute/compute_fluxtube_growth_rate.m +++ b/matlab/compute/compute_fluxtube_growth_rate.m @@ -1,72 +1,73 @@ function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT) % Remove AA part if DATA.Nx > 1 ikxnz = abs(DATA.PHI(:,1,1,1)) > 0; else ikxnz = abs(DATA.PHI(:,1,1,1)) > -1; end ikynz = abs(DATA.PHI(1,:,1,1)) > 0; phi = DATA.PHI(ikxnz,ikynz,:,:); t = DATA.Ts3D; [~,its] = min(abs(t-TRANGE(1))); [~,ite] = min(abs(t-TRANGE(end))); w_ky = zeros(sum(ikynz),ite-its); ce = zeros(sum(ikynz),ite-its); is = 1; for it = its+1:ite phi_n = phi(:,:,:,it); phi_nm1 = phi(:,:,:,it-1); dt = t(it)-t(it-1); ZS = sum(sum(phi_nm1,1),3); wl = log(phi_n./phi_nm1)/dt; w_ky(:,is) = sum(sum(wl.*phi_nm1,1),3)./ZS; for iky = 1:numel(w_ky(:,is)) ce(iky,is) = abs(sum(sum(abs(w_ky(iky,is)-wl(:,iky,:)).^2.*phi_nm1(:,iky,:),1),3)./ZS(:,iky,:)); end is = is + 1; end [kys, Is] = sort(DATA.ky(ikynz)); linear_gr.trange = t(its:ite); linear_gr.g_ky = real(w_ky(Is,:)); linear_gr.w_ky = imag(w_ky(Is,:)); linear_gr.ce = abs(ce); linear_gr.ky = kys; if PLOT >0 figure plot(linear_gr.ky,linear_gr.g_ky(:,end),'DisplayName','$Re(\omega_{k_y})$'); hold on; plot(linear_gr.ky,linear_gr.w_ky(:,end),'DisplayName','$Im(\omega_{k_y})$'); hold on; plot(linear_gr.ky,linear_gr.ce (:,end),'DisplayName','$\epsilon$'); hold on; xlim([min(linear_gr.ky) max(linear_gr.ky)]); xlabel('$k_y$'); legend('show'); + title(DATA.param_title); if PLOT > 1 [YY,XX] = meshgrid(kys,t(its+1:ite)); figure subplot(311) % imagesc(t(its:ite),kys,real(w_ky)); set(gca,'YDir','normal'); pclr = pcolor(XX',YY',real(w_ky)); set(pclr, 'edgecolor','none'); xlabel('$t$'); ylabel('$k_y$'); title('Re($\omega$)') subplot(312) pclr = pcolor(XX',YY',imag(w_ky)); set(pclr, 'edgecolor','none'); xlabel('$t$'); ylabel('$k_y$'); title('Im($\omega$)') subplot(313) pclr = pcolor(XX',YY',abs(w_ky)); set(pclr, 'edgecolor','none'); xlabel('$t$'); ylabel('$k_y$'); title('|$\omega$|') end end end \ No newline at end of file diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m index 95af8e7..5dc0088 100644 --- a/matlab/compute/mode_growth_meter.m +++ b/matlab/compute/mode_growth_meter.m @@ -1,99 +1,99 @@ function [FIGURE] = mode_growth_meter(DATA,OPTIONS) NORMALIZED = OPTIONS.NORMALIZED; Nma = OPTIONS.NMA; %Number moving average t = OPTIONS.TIME; -iz = 1; +iz = OPTIONS.iz; [~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2))); FRAMES = zeros(size(OPTIONS.TIME)); for i = 1:numel(OPTIONS.TIME) [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D)); end FIGURE.fig = figure; set(gcf, 'Position', [100 100 1200 700]) FIGURE.FIGNAME = 'mode_growth_meter'; for i = 1:2 MODES_SELECTOR = i; %(1:Zonal, 2: NZonal, 3: ky=kx) if MODES_SELECTOR == 1 if NORMALIZED plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma); else plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma); end kstr = 'k_x'; k = DATA.kx; MODESTR = 'Zonal modes'; elseif MODES_SELECTOR == 2 if NORMALIZED plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma); else plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma); end kstr = 'k_y'; k = DATA.ky; MODESTR = 'NZ modes'; elseif MODES_SELECTOR == 3 if NORMALIZED plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),Nma); else plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),Nma); end kstr = 'k_y=k_x'; k = DATA.ky; MODESTR = 'Diag modes'; end MODES = 1:numel(k); % MODES = zeros(size(OPTIONS.K2PLOT)); % for i = 1:numel(OPTIONS.K2PLOT) % [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k)); % end % plt = @(x,ik) abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))); gamma = MODES; amp = MODES; i_=1; for ik = MODES gr = polyfit(t',log(plt(DATA.PHI,ik)),1); gamma(i_) = gr(1); amp(i_) = gr(2); i_=i_+1; end %plot subplot(2,3,1+3*(i-1)) [YY,XX] = meshgrid(t,fftshift(k,1)); pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on; set(gca,'YDir','normal') % xlim([t(1) t(end)]); %ylim([1e-5 1]) xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$'); title(MODESTR) subplot(2,3,2+3*(i-1)) mod2plot = [2:OPTIONS.NMODES+1]; for i_ = mod2plot semilogy(t,plt(DATA.PHI,MODES(i_))); hold on; % semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k') end if MODES_SELECTOR == 1 semilogy(t,plt(DATA.PHI,ikzf),'--k'); end xlim([t(1) t(end)]); %ylim([1e-5 1]) xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']); title('Measure time window') subplot(2,3,3+3*(i-1)) plot(k(MODES),gamma); hold on; plot(k(MODES(mod2plot)),gamma(mod2plot),'x') if MODES_SELECTOR == 1 plot(k(ikzf),gamma(ikzf),'ok'); end xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$'); title('Growth rates') end end \ No newline at end of file diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m index 3eddddf..6d2dfa7 100644 --- a/matlab/load/quick_gene_load.m +++ b/matlab/load/quick_gene_load.m @@ -1,46 +1,46 @@ %% gyroLES plot % genepath = '/misc/gene_results'; % % simpath = '/NL_Zpinch_Kn_1.7_eta_0.25_nuSG_1e-1_gyroLES/'; % simpath = '/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_gyroLES/'; % filename = 'GyroLES_ions.dat'; % toload = [genepath simpath filename]; % data = readtable(toload); % t = data.Var1; % mu_x = data.Var2; % mu_y = data.Var3; % figure; % plot(t,mu_x); hold on % plot(t,mu_y); % legend('mu_x','mu_y'); xlabel('time'); ylabel('Hyper-diff'); % title('gyroLES results for H&P fig. 2c') %% Plot linear results path = '/home/ahoffman/gene/linear_zpinch_results/'; % fname ='tmp.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.047_32x16.txt'; % fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.047_32x16.txt'; % fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0235_32x16.txt'; % fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_32x16.txt'; % fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_64x32.txt'; % fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt'; % fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt'; % fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSGDK_0.0047_64x32.txt'; -fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt'; +% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt'; -% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt'; +fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt'; % fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.235_64x32.txt'; % fname ='GENE_LIN_Kn_1.7_KT_0.425_nuSG_0.235_64x32.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.475_nu_0_mu_5e-2.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_64x32.txt'; % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_32x16.txt'; % fname ='GENE_LIN_Kn_2.0_KT_0.5_nu_0_32x16.txt'; % fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt'; % fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt'; % fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt'; data_ = load([path,fname]); figure plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE'); \ No newline at end of file diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m index ed66da6..4c7ee67 100644 --- a/matlab/plot/create_film.m +++ b/matlab/plot/create_film.m @@ -1,116 +1,120 @@ function create_film(DATA,OPTIONS,format) %% Plot options FPS = 30; DELAY = 1/FPS; BWR = 1; NORMALIZED = 1; if ~strcmp(OPTIONS.PLAN,'sx') T = DATA.Ts3D; else T = DATA.Ts5D; end %% Processing -toplot = process_field(DATA,OPTIONS); - +switch OPTIONS.PLAN + case 'RZ' + toplot = poloidal_plot(DATA,OPTIONS); + otherwise + toplot = process_field(DATA,OPTIONS); +end %% FILENAME = [DATA.localdir,toplot.FILENAME,format]; XNAME = ['$',toplot.XNAME,'$']; YNAME = ['$',toplot.YNAME,'$']; FIELDNAME = ['$',toplot.FIELDNAME,'$']; FIELD = toplot.FIELD; X = toplot.X; Y = toplot.Y; FRAMES = toplot.FRAMES; switch format case '.avi' vidfile = VideoWriter(FILENAME,'Uncompressed AVI'); vidfile.FrameRate = FPS; open(vidfile); end % Set colormap boundaries -hmax = max(max(max(FIELD(:,:,FRAMES)))); -hmin = min(min(min(FIELD(:,:,FRAMES)))); +hmax = max(max(max(FIELD(:,:,:)))); +hmin = min(min(min(FIELD(:,:,:)))); scale = -1; flag = 0; if hmax == hmin disp('Warning : h = hmin = hmax = const') else % Setup figure frame fig = figure('Color','white','Position', toplot.DIMENSIONS.*[1 1 1.2 1]); if ~strcmp(OPTIONS.PLAN,'sx') pcolor(X,Y,FIELD(:,:,1)); % to set up if BWR colormap(bluewhitered) end axis tight manual % this ensures that getframe() returns a consistent size if toplot.INTERP shading interp; end else contour(toplot.X,toplot.Y,FIELD(:,:,1),128); end in = 1; nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); - for n = FRAMES % loop over selected frames + for n = 1:numel(FRAMES) % loop over selected frames scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize if ~strcmp(OPTIONS.PLAN,'sx') if NORMALIZED pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\ caxis([-1,1]); else pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\ if CONST_CMAP caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map else caxis([-1,1]*scale); % adaptive color map end title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... ,', scale = ',sprintf('%.1e',scale)]); end if toplot.INTERP shading interp; end set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT); if BWR colormap(bluewhitered) end else % show velocity distr. contour(toplot.X,toplot.Y,FIELD(:,:,n)/scale,128); end - title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... + title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(FRAMES(n))))... ,', scaling = ',sprintf('%.1e',scale)]); xlabel(XNAME); ylabel(YNAME); %colorbar; drawnow % Capture the plot as an image frame = getframe(fig); switch format case '.gif' im = frame2im(frame); [imind,cm] = rgb2ind(im,32); % Write to the GIF File if in == 1 imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf); else imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY); end case '.avi' writeVideo(vidfile,frame); otherwise disp('Unknown format'); break 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(' ') switch format case '.gif' disp(['Gif saved @ : ',FILENAME]) case '.avi' disp(['Video saved @ : ',FILENAME]) close(vidfile); end end end diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m index 4a780f3..278fe02 100644 --- a/matlab/plot/photomaton.m +++ b/matlab/plot/photomaton.m @@ -1,46 +1,51 @@ function [ FIGURE ] = photomaton( DATA,OPTIONS ) %UNTITLED5 Summary of this function goes here %% Processing -toplot = process_field(DATA,OPTIONS); +switch OPTIONS.PLAN + case 'RZ' + toplot = poloidal_plot(DATA,OPTIONS); + otherwise + toplot = process_field(DATA,OPTIONS); +end FNAME = toplot.FILENAME; FRAMES = toplot.FRAMES; Nframes= numel(FRAMES); Nrows = ceil(Nframes/4); Ncols = ceil(Nframes/Nrows); % TNAME = []; FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows]) for i_ = 1:numel(FRAMES) subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))]; if ~strcmp(OPTIONS.PLAN,'kxky') - scale = max(max(abs(toplot.FIELD(:,:,FRAMES(i_))))); % Scaling to normalize + scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize else scale = 1; end if ~strcmp(OPTIONS.PLAN,'sx') tshot = DATA.Ts3D(FRAMES(i_)); - pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none'); + pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none'); if OPTIONS.AXISEQUAL pbaspect(toplot.ASPECT) end if ~strcmp(OPTIONS.PLAN,'kxky') caxis([-1,1]); colormap(bluewhitered); end if OPTIONS.INTERP shading interp; end else - contour(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale,128); + contour(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale,128); % pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none'); shading interp tshot = DATA.Ts5D(FRAMES(i_)); end xlabel(toplot.XNAME); ylabel(toplot.YNAME); % if i_ > 1; set(gca,'ytick',[]); end; title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',scale)]); end legend(['$',toplot.FIELDNAME,'$']); FIGURE.FIGNAME = [FNAME,'_snaps',TNAME]; end diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m new file mode 100644 index 0000000..f209aea --- /dev/null +++ b/matlab/plot/plot_ballooning.m @@ -0,0 +1,53 @@ +function [FIG] = plot_ballooning(data,options) + + [~,it] = min(abs(data.Ts3D - options.time_2_plot)); + phi_real=(real(data.PHI(:,:,:,it))); + phi_imag=(imag(data.PHI(:,:,:,it))); + % Apply baollooning tranform + for iky=2 + dims = size(phi_real); + phib_real = zeros( dims(1)*dims(3) ,1); + phib_imag= phib_real; + b_angle = phib_real; + + midpoint = floor((dims(1)*dims(3) )/2)+1; + for ip =1: dims(1) + start_ = (ip -1)*dims(3) +1; + end_ = ip*dims(3); + phib_real(start_:end_) = phi_real(ip,iky,:); + phib_imag(start_:end_) = phi_imag(ip,iky, :); + end + + % Define ballooning angle + Nkx = numel(data.kx)-1; coordz = data.z; + idx = -Nkx:1:Nkx; + for ip =1: dims(1) + for iz=1:dims(3) + ii = dims(3)*(ip -1) + iz; + b_angle(ii) =coordz(iz) + 2*pi*idx(ip); + end + end + + phib = phib_real(:) + 1i * phib_imag(:); + % normalize real and imaginary parts at chi =0 + if options.normalized + [~,idxLFS] = min(abs(b_angle -0)); + normalization = phib( idxLFS); + else + normalization = 1; + end + phib_norm = phib / normalization ; + phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS); + phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS); + + + FIG.fig = figure; hold all; + plot(b_angle/pi, phib_real_norm,'-b'); + plot(b_angle/pi, phib_imag_norm ,'-r'); + plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k'); + legend('real','imag','norm') + xlabel('$\chi / \pi$') + ylabel('$\phi_B (\chi)$'); + title(['HeLaZ ballooning, t=',num2str(data.Ts3D(it))]); + end +end diff --git a/matlab/plot/poloidal_plot.m b/matlab/plot/poloidal_plot.m new file mode 100644 index 0000000..66775a9 --- /dev/null +++ b/matlab/plot/poloidal_plot.m @@ -0,0 +1,73 @@ +function [TOPLOT] = poloidal_plot(DATA,OPTIONS) +sq = @(x) squeeze(x); +%% Geometry +r0 = DATA.a; +q0 = DATA.Q0; +Cy = r0/q0; +Ly = max(DATA.y) - min(DATA.y); +n0 = Cy*2*pi/Ly; +z_ = DATA.z; +%% grid sizes +nkx = DATA.Nkx; nky = DATA.Nky; nz = DATA.Nz; +nx = DATA.Nx; ny = DATA.Ny; +%% Time and frames +FRAMES = zeros(size(OPTIONS.TIME)); +for i = 1:numel(OPTIONS.TIME) + [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D)); +end +FRAMES = unique(FRAMES); + +%% Field to plot +switch OPTIONS.NAME + case '\phi' + FIELD = DATA.PHI; + NAME = 'phi'; + case 'n_i' + FIELD = DATA.DENS_I; + NAME = 'ni'; +end + +%% Build poloidal plot (tor angle = 0) +FRZ = zeros(DATA.Nx,DATA.Nz, numel(FRAMES)); +w_ = zeros(nky,nz); +for iz = 1:nz + for iky = 1:ny + w_(iky,iz) = exp(1i*iky*n0*(q0*z_(iz))); + end +end + + +for it = 1:numel(FRAMES) + field_ = squeeze(FIELD(:,:,:,FRAMES(it))); + spectrumKxKyZ = zeros(nkx,nky,nz); + spectrumKxKyZ(1:nkx,:,:) = field_; + spectrumKxKyZ((nkx+1):nx,1,:) = conj(field_(nkx:-1:2,1,:)); + spectrumKxKyZ((nkx+1):nx,2:ny,:) = conj(field_(nkx:-1:2,ny:-1:2,:)); + + Fxkyz = ny*ifft(spectrumKxKyZ,[],1); + for ix = 1:DATA.Nx + FRZ(ix,:,it) = sq(sum(sq(Fxkyz(ix,:,:)).*w_(:,:),1)); + end +end +FRZ = real(FRZ); +%grid +[Y,X] = meshgrid(DATA.z,DATA.x); +X__ = (X+DATA.a).*cos(Y); +Y__ = (X+DATA.a).*sin(Y); +X = X__; +Y = Y__; +%% output + +TOPLOT.FIELD = FRZ; +TOPLOT.FRAMES = FRAMES; +TOPLOT.X = (X); +TOPLOT.Y = (Y); +TOPLOT.FIELDNAME = NAME; +TOPLOT.XNAME = '$R$'; +TOPLOT.YNAME = '$Z$'; +TOPLOT.FILENAME = [NAME,'_poloidal_phieq0']; +TOPLOT.DIMENSIONS= [100, 100, 500, 500]; +TOPLOT.ASPECT = [1 1 1]; +TOPLOT.INTERP = OPTIONS.INTERP; + +end \ No newline at end of file diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m index 5dc78fa..6e80a71 100644 --- a/matlab/plot/show_moments_spectrum.m +++ b/matlab/plot/show_moments_spectrum.m @@ -1,142 +1,157 @@ function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS ) Pi = DATA.Pi; Ji = DATA.Ji; Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5); Nipj = squeeze(Nipj); +if DATA.K_E Pe = DATA.Pe; Je = DATA.Je; Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5); Nepj = squeeze(Nepj); +end FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS]; set(gcf, 'Position', [100 50 1000 400]) if ~OPTIONS.ST t0 = OPTIONS.TIME(1); t1 = OPTIONS.TIME(end); [~,it0] = min(abs(t0-DATA.Ts5D)); [~,it1] = min(abs(t1-DATA.Ts5D)); Nipj = mean(Nipj(:,:,it0:it1),3); + if DATA.K_E Nepj = mean(Nepj(:,:,it0:it1),3); + end if numel(OPTIONS.TIME) == 1 TITLE=['$t=',num2str(OPTIONS.TIME),'$']; else TITLE=['$t\in[',num2str(t0),',',num2str(t1),']$']; end Nipj = squeeze(Nipj); - Nepj = squeeze(Nepj); ymini = min(Nipj); ymaxi = max(Nipj); + + if DATA.K_E + Nepj = squeeze(Nepj); ymine = min(Nepj); ymaxe = max(Nepj); - ymin = min([ymini ymine]); ymax = max([ymaxi ymaxe]); - - + ymin = min([ymini ymine]); + end + if DATA.K_E subplot(121) + end if ~OPTIONS.P2J for ij = 1:numel(Ji) name = ['$j=',num2str(Ji(ij)),'$']; semilogy(Pi,Nipj(:,ij),'o-','DisplayName',name); hold on; end xlabel('$p$'); else for ij = 1:numel(Ji) name = ['$j=',num2str(Ji(ij)),'$']; semilogy(Pi+2*Ji(ij),Nipj(:,ij),'o-','DisplayName',name); hold on; end xlabel('$p+2j$') end ylabel(['$\sum_{kx,ky}|N_i^{pj}|$']); legend('show'); title([TITLE,' He-La ion spectrum']); - + + if DATA.K_E subplot(122) if ~OPTIONS.P2J for ij = 1:numel(Je) name = ['$j=',num2str(Je(ij)),'$']; semilogy(Pe,Nepj(:,ij),'o-','DisplayName',name); hold on; end xlabel('p'); else for ij = 1:numel(Je) name = ['$j=',num2str(Je(ij)),'$']; semilogy(Pe+2*Je(ij),Nepj(:,ij),'o-','DisplayName',name); hold on; end xlabel('$p+2j$') end ylabel(['$\sum_{kx,ky}|N_e^{pj}|$']); legend('show'); title([TITLE,' He-La elec. spectrum']); + end else [JJ,PP] = meshgrid(Ji,Pi); P2Ji = PP + 2*JJ; % form an axis of velocity ordered moments p2ji = unique(P2Ji); % weights to average weights = zeros(size(p2ji)); % space time of moments amplitude wrt to p+2j degree Ni_ST = zeros(numel(p2ji),numel(DATA.Ts5D)); % fill the st diagramm by averaging every p+2j=n moments together for ip = 1:numel(Pi) for ij = 1:numel(Ji) [~,ip2j] = min(abs(p2ji-(Pi(ip)+2*Ji(ij)))); Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:))); weights(ip2j) = weights(ip2j) + 1; end end % doing the average for ip2j = 1:numel(p2ji) Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j); end + if DATA.K_E % same for electrons!! [JJ,PP] = meshgrid(Je,Pe); P2Je = PP + 2*JJ; % form an axis of velocity ordered moments p2je = unique(P2Je); % weights to average weights = zeros(size(p2je)); % space time of moments amplitude wrt to p+2j degree Ne_ST = zeros(numel(p2je),numel(DATA.Ts5D)); % fill the st diagramm by averaging every p+2j=n moments together for ip = 1:numel(Pe) for ij = 1:numel(Je) [~,ip2j] = min(abs(p2ji-(Pe(ip)+2*Je(ij)))); Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:))); weights(ip2j) = weights(ip2j) + 1; end end % doing the average for ip2j = 1:numel(p2ji) Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j); end + end % plots % scaling % plt = @(x,ip2j) x./max(max(x)); if OPTIONS.NORMALIZED plt = @(x,ip2j) x(ip2j,:)./max(x(ip2j,:)); else plt = @(x,ip2j) x; end + if DATA.K_E subplot(2,1,1) + end imagesc(DATA.Ts5D,p2ji,plt(Ni_ST,1:numel(p2ji))); set(gca,'YDir','normal') % pclr = pcolor(XX,YY,plt(Ni_ST)); % set(pclr, 'edgecolor','none'); hold on; xlabel('$t$'); ylabel('$p+2j$') title('$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$') + if DATA.K_E subplot(2,1,2) imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ji))); set(gca,'YDir','normal') % pclr = pcolor(XX,YY,plt(Ne_ST)); % set(pclr, 'edgecolor','none'); hold on; xlabel('$t$'); ylabel('$p+2j$') title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$') suptitle(DATA.param_title); + end end end diff --git a/matlab/process_field.m b/matlab/process_field.m index f19f1b2..5762fc8 100644 --- a/matlab/process_field.m +++ b/matlab/process_field.m @@ -1,450 +1,451 @@ function [ TOPLOT ] = process_field( DATA, OPTIONS ) %UNTITLED4 Summary of this function goes here % Detailed explanation goes here %% Setup time %% FRAMES = zeros(size(OPTIONS.TIME)); for i = 1:numel(OPTIONS.TIME) [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D)); end + %% Setup the plot geometry [KY,~] = meshgrid(DATA.ky,DATA.kx); directions = {'x','y','z'}; -Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(OPTIONS.TIME); +Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES); POLARPLOT = OPTIONS.POLARPLOT; LTXNAME = OPTIONS.NAME; switch OPTIONS.PLAN case 'xy' XNAME = '$x$'; YNAME = '$y$'; [Y,X] = meshgrid(DATA.y,DATA.x); REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; case 'xz' XNAME = '$x$'; YNAME = '$z$'; [Y,X] = meshgrid(DATA.z,DATA.x); REALP = 1; COMPDIM = 2; SCALE = 0; case 'yz' XNAME = '$y$'; YNAME = '$z$'; [Y,X] = meshgrid(DATA.z,DATA.y); REALP = 1; COMPDIM = 1; SCALE = 0; case 'kxky' XNAME = '$k_x$'; YNAME = '$k_y$'; [Y,X] = meshgrid(DATA.ky,DATA.kx); REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; case 'kxz' XNAME = '$k_x$'; YNAME = '$z$'; [Y,X] = meshgrid(DATA.z,DATA.kx); REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0; case 'kyz' XNAME = '$k_y$'; YNAME = '$z$'; [Y,X] = meshgrid(DATA.z,DATA.ky); REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0; case 'sx' XNAME = '$v_\parallel$'; YNAME = '$\mu$'; [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR); REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0; end Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y))); Lmin_ = min([Xmax_,Ymax_]); Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2; Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2; DIMENSIONS = [100, 100, 400*Rx, 400*Ry]; ASPECT = [Rx, Ry, 1]; % Polar grid POLARNAME = []; if POLARPLOT POLARNAME = 'polar'; X__ = (X+DATA.a).*cos(Y); Y__ = (X+DATA.a).*sin(Y); X = X__; Y = Y__; XNAME='X'; YNAME='Z'; DIMENSIONS = [100, 100, 500, 500]; ASPECT = [1,1,1]; sz = size(X); FIELD = zeros(sz(1),sz(2),Nt); else sz = size(X); FIELD = zeros(sz(1),sz(2),Nt); end %% Process the field to plot switch OPTIONS.COMP case 'avg' compr = @(x) mean(x,COMPDIM); COMPNAME = ['avg',directions{COMPDIM}]; FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}]; case 'max' compr = @(x) max(x,[],COMPDIM); COMPNAME = ['max',directions{COMPDIM}]; FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME]; otherwise switch COMPDIM case 1 i = OPTIONS.COMP; compr = @(x) x(i,:,:); if REALP COMPNAME = sprintf(['x=','%2.1f'],DATA.x(i)); else COMPNAME = sprintf(['k_x=','%2.1f'],DATA.kx(i)); end FIELDNAME = [LTXNAME,'(',COMPNAME,')']; case 2 i = OPTIONS.COMP; compr = @(x) x(:,i,:); if REALP COMPNAME = sprintf(['y=','%2.1f'],DATA.y(i)); else COMPNAME = sprintf(['k_y=','%2.1f'],DATA.ky(i)); end FIELDNAME = [LTXNAME,'(',COMPNAME,')']; case 3 i = OPTIONS.COMP; compr = @(x) x(:,:,i); COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i)); FIELDNAME = [LTXNAME,'(',COMPNAME,')']; end end switch REALP case 1 % Real space plot INTERP = OPTIONS.INTERP; process = @(x) real(fftshift(ifft2(x,Nx,Ny))); shift_x = @(x) x; shift_y = @(x) x; case 0 % Frequencies plot INTERP = 0; switch COMPDIM case 1 process = @(x) abs(fftshift(x,2)); shift_x = @(x) fftshift(x,1); shift_y = @(x) fftshift(x,1); case 2 process = @(x) abs((x)); shift_x = @(x) (x); shift_y = @(x) (x); case 3 process = @(x) abs(fftshift(x,2)); shift_x = @(x) fftshift(x,2); shift_y = @(x) fftshift(x,2); end end switch OPTIONS.NAME case '\phi' NAME = 'phi'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.PHI(:,:,:,it))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it)))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'N_i^{00}' NAME = 'Ni00'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.Ni00(:,:,:,it))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.Ni00(:,:,:,FRAMES(it)))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'n_e' NAME = 'ne'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.DENS_E(:,:,:,it))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it)))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'k^2n_e' NAME = 'k2ne'; [KY, KX] = meshgrid(DATA.ky, DATA.kx); if COMPDIM == 3 - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:DATA.Nz - tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it))); + tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'n_e^{NZ}' NAME = 'neNZ'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.DENS_E(:,:,:,it).*(KY~=0))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it)).*(KY~=0))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it).*(KY~=0))); + tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it)).*(KY~=0))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'n_i' NAME = 'ni'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.DENS_I(:,:,:,it))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it)))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'T_i' NAME = 'Ti'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.TEMP_I(:,:,:,it))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.TEMP_I(:,:,:,FRAMES(it)))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,it))); + tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,FRAMES(it)))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'n_i^{NZ}' NAME = 'niNZ'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.DENS_I(:,:,:,it).*(KY~=0))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it)).*(KY~=0))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it).*(KY~=0))); + tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it)).*(KY~=0))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case '\phi^{NZ}' NAME = 'phiNZ'; if COMPDIM == 3 - for it = FRAMES - tmp = squeeze(compr(DATA.PHI(:,:,:,it).*(KY~=0))); + for it = 1:numel(FRAMES) + tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it)).*(KY~=0))); FIELD(:,:,it) = squeeze(process(tmp)); end else if REALP tmp = zeros(Nx,Ny,Nz); else tmp = zeros(DATA.Nkx,DATA.Nky,Nz); end - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it).*(KY~=0))); + tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it)).*(KY~=0))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'v_y' NAME = 'vy'; [~, KX] = meshgrid(DATA.ky, DATA.kx); vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES)); - for it = FRAMES % Compute the 3D real transport for each frame + for it = 1:numel(FRAMES) for iz = 1:DATA.Nz - vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)); + vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny)); end end if REALP % plot in real space - for it = FRAMES + for it = 1:numel(FRAMES) FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it))); end else % Plot spectrum process = @(x) abs(fftshift(x,2)); tmp = zeros(DATA.Nkx,DATA.Nky,Nz); - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,it)))); + tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 'v_x' NAME = 'vx'; [KY, ~] = meshgrid(DATA.ky, DATA.kx); vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES)); - for it = FRAMES % Compute the 3D real transport for each frame + for it = 1:numel(FRAMES) for iz = 1:DATA.Nz - vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)); + vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny)); end end if REALP % plot in real space - for it = FRAMES + for it = 1:numel(FRAMES) FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it))); end else % Plot spectrum process = @(x) abs(fftshift(x,2)); tmp = zeros(DATA.Nkx,DATA.Nky,Nz); - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,it)))); + tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case 's_y' NAME = 'sy'; [~, KX] = meshgrid(DATA.ky, DATA.kx); vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES)); - for it = FRAMES % Compute the 3D real transport for each frame + for it = 1:numel(FRAMES) for iz = 1:DATA.Nz - vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)); + vE(:,:,iz,it) = real(ifft2(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny)); end end if REALP % plot in real space - for it = FRAMES + for it = 1:numel(FRAMES) FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it))); end else % Plot spectrum process = @(x) abs(fftshift(x,2)); tmp = zeros(DATA.Nkx,DATA.Nky,Nz); - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) - tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,it)))); + tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it))))); end FIELD(:,:,it) = squeeze(compr(tmp)); end end case '\Gamma_x' NAME = 'Gx'; [KY, ~] = meshgrid(DATA.ky, DATA.kx); Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES)); - for it = FRAMES % Compute the 3D real transport for each frame + for it = 1:numel(FRAMES) for iz = 1:DATA.Nz - Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))... - .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny))); + Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Nx,DATA.Ny)))... + .*real((ifft2(DATA.DENS_I(:,:,iz,FRAMES(it)),DATA.Nx,DATA.Ny))); end end if REALP % plot in real space - for it = FRAMES + for it = 1:numel(FRAMES) FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it))); end else % Plot spectrum process = @(x) abs(fftshift(x,2)); shift_x = @(x) fftshift(x,2); shift_y = @(x) fftshift(x,2); tmp = zeros(DATA.Nx,DATA.Ny,Nz); - for it = FRAMES + for it = 1:numel(FRAMES) for iz = 1:numel(DATA.z) tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny))); end FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:))); end end case 'f_i' shift_x = @(x) x; shift_y = @(y) y; NAME = 'fi'; OPTIONS.SPECIE = 'i'; [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D)); [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D)); dit_ = 1;%ceil((it1_-it0_+1)/10); FRAMES = it0_:dit_:it1_; iz = 1; - for it = FRAMES - OPTIONS.T = DATA.Ts5D(it); + for it = 1:numel(FRAMES) + OPTIONS.T = DATA.Ts5D(FRAMES(it)); [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS); end case 'f_e' shift_x = @(x) x; shift_y = @(y) y; NAME = 'fe'; OPTIONS.SPECIE = 'e'; [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D)); [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D)); dit_ = 1;%ceil((it1_-it0_+1)/10); FRAMES = it0_:dit_:it1_; iz = 1; - for it = FRAMES - OPTIONS.T = DATA.Ts5D(it); + for it = 1:numel(FRAMES) + OPTIONS.T = DATA.Ts5D(FRAMES(it)); [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS); end end TOPLOT.FIELD = FIELD; TOPLOT.FRAMES = FRAMES; TOPLOT.X = shift_x(X); TOPLOT.Y = shift_y(Y); TOPLOT.FIELDNAME = FIELDNAME; TOPLOT.XNAME = XNAME; TOPLOT.YNAME = YNAME; TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME]; TOPLOT.DIMENSIONS= DIMENSIONS; TOPLOT.ASPECT = ASPECT; TOPLOT.FRAMES = FRAMES; TOPLOT.INTERP = INTERP; end diff --git a/matlab/profiler.m b/matlab/profiler.m index 2988e22..a8184e1 100644 --- a/matlab/profiler.m +++ b/matlab/profiler.m @@ -1,96 +1,103 @@ %% load profiling % filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00); -% filename_ = ['/misc/HeLaZ_outputs',filename(3:end)]; -filename_ = ['/home/ahoffman/HeLaZ',filename(3:end)]; -CPUTIME = double(h5readatt(filename_,'/data/input','cpu_time')); -DT_SIM = h5readatt(filename_,'/data/input','dt'); +% outfilename = ['/misc/HeLaZ_outputs',filename(3:end)]; +outfilename = data.outfilenames{end}; +CPUTIME = double(h5readatt(outfilename,'/data/input','cpu_time')); +DT_SIM = h5readatt(outfilename,'/data/input','dt'); -rhs_Tc = h5read(filename_,'/profiler/Tc_rhs'); -adv_field_Tc = h5read(filename_,'/profiler/Tc_adv_field'); -ghost_Tc = h5read(filename_,'/profiler/Tc_ghost'); -clos_Tc = h5read(filename_,'/profiler/Tc_clos'); -coll_Tc = h5read(filename_,'/profiler/Tc_coll'); -poisson_Tc = h5read(filename_,'/profiler/Tc_poisson'); -Sapj_Tc = h5read(filename_,'/profiler/Tc_Sapj'); -checkfield_Tc= h5read(filename_,'/profiler/Tc_checkfield'); -diag_Tc = h5read(filename_,'/profiler/Tc_diag'); -step_Tc = h5read(filename_,'/profiler/Tc_step'); -Ts0D = h5read(filename_,'/profiler/time'); +rhs_Tc = h5read(outfilename,'/profiler/Tc_rhs'); +adv_field_Tc = h5read(outfilename,'/profiler/Tc_adv_field'); +ghost_Tc = h5read(outfilename,'/profiler/Tc_ghost'); +clos_Tc = h5read(outfilename,'/profiler/Tc_clos'); +coll_Tc = h5read(outfilename,'/profiler/Tc_coll'); +poisson_Tc = h5read(outfilename,'/profiler/Tc_poisson'); +Sapj_Tc = h5read(outfilename,'/profiler/Tc_Sapj'); +checkfield_Tc= h5read(outfilename,'/profiler/Tc_checkfield'); +diag_Tc = h5read(outfilename,'/profiler/Tc_diag'); +process_Tc = h5read(outfilename,'/profiler/Tc_process'); +step_Tc = h5read(outfilename,'/profiler/Tc_step'); +Ts0D = h5read(outfilename,'/profiler/time'); + +N_T = 11; missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ... - -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc; + -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc-process_Tc; total_Tc = step_Tc; TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(ghost_Tc);... diff(clos_Tc); diff(coll_Tc); diff(poisson_Tc); diff(Sapj_Tc); ... - diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)]; -TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/10,10]); + diff(checkfield_Tc); diff(diag_Tc); diff(process_Tc); diff(missing_Tc)]; +TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]); TIME_PER_STEP = sum(TIME_PER_FCT,2); TIME_PER_CPU = trapz(Ts0D(2:end),TIME_PER_STEP); rhs_Ta = mean(diff(rhs_Tc)); adv_field_Ta = mean(diff(adv_field_Tc)); ghost_Ta = mean(diff(ghost_Tc)); clos_Ta = mean(diff(clos_Tc)); coll_Ta = mean(diff(coll_Tc)); poisson_Ta = mean(diff(poisson_Tc)); Sapj_Ta = mean(diff(Sapj_Tc)); checkfield_Ta = mean(diff(checkfield_Tc)); +process_Ta = mean(diff(process_Tc)); diag_Ta = mean(diff(diag_Tc)); NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM; %% Plots if 1 %% Area plot fig = figure; - +% colors = rand(N_T,3); +colors = lines(N_T); p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none'); -for i = 1:10; p1(i).FaceColor = rand(1,3); end; -legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing') +for i = 1:N_T; p1(i).FaceColor = colors(i,:); end; +legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing') xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]') xlim([Ts0D(2),Ts0D(end)]); title(sprintf('Proc. 1, total sim. time ~%.0f [h]',CPUTIME/3600)) hold on FIGNAME = 'profiler'; -save_figure +% save_figure else %% Normalized Area plot fig = figure; p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none', 'FaceColor','flat'); -for i = 1:10; p1(i).FaceColor = rand(1,3); end; +for i = 1:N_T; p1(i).FaceColor = rand(1,3); end; legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing') xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]') ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]); hold on yyaxis right p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0); ylabel('Step Comp. Time [s]') ylim([0,1.1*max(diff(total_Tc))]) set(gca,'ycolor','r') FIGNAME = 'profiler'; -save_figure +% save_figure end if 0 %% Histograms fig = figure; histogram(diff(rhs_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(adv_field_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(ghost_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(clos_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(coll_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(poisson_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(Sapj_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on +histogram(diff(process_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(checkfield_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on histogram(diff(diag_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on grid on; -legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing') +legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Process','Check+sym', 'Diagnos.', 'Missing') xlabel('Step Comp. Time [s]'); ylabel('') +set(gca,'Xscale','log') FIGNAME = 'profiler'; -save_figure +% save_figure end \ No newline at end of file diff --git a/matlab/setup.m b/matlab/setup.m index 00269b3..5233fad 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,174 +1,175 @@ -%% ________________________________________________________________________ +%% _______________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nx = NX; % x grid resolution GRID.Lx = LX; % x length GRID.Ny = NY; % y '' GRID.Ly = LY; % y '' GRID.Nz = NZ; % z resolution +GRID.Npol = NPOL; % z resolution if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end; % Geometry GEOM.geom = ['''',GEOMETRY,'''']; GEOM.q0 = Q0; % q factor GEOM.shear = SHEAR; % shear GEOM.eps = EPS; % inverse aspect ratio % Model parameters MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; MODEL.LINEARITY = ['''',LINEARITY,'''']; MODEL.KIN_E = KIN_E; if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end; MODEL.mu_x = MU_X; MODEL.mu_y = MU_Y; MODEL.mu_z = 0; 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 = SIGMA_E; 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.K_n = K_N; % source term kappa for HW MODEL.K_T = K_T; % Temperature MODEL.K_E = K_E; % Electric MODEL.GradB = GRADB; % Magnetic gradient MODEL.CurvB = CURVB; % Magnetic curvature MODEL.lambdaD = LAMBDAD; % Collision parameters COLL.collision_model = ['''',CO,'''']; if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end; if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end; COLL.mat_file = '''null'''; switch CO case 'SG' COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5'''; case 'LR' COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''; case 'LD' % COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5'''; COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''; end % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; INITIAL.INIT_OPT = ['''',INIT_OPT,'''']; INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,'''']; INITIAL.init_background = BCKGD0; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; % Naming and creating input file CONAME = CO; if GKCO CONAME = [CONAME,'GK']; else CONAME = [CONAME,'DK']; end if ~ABCO CONAME = [CONAME,'aa']; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end % Hermite-Laguerre degrees naming if (PMAXE == PMAXI) && (JMAXE == JMAXI) HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)]; else HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),... '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)]; end % temp. dens. drives drives_ = []; if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end; if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end; % collision coll_ = ['_nu_%0.0e_',CONAME]; coll_ = sprintf(coll_,NU); % nonlinear lin_ = []; if ~LINEARITY; lin_ = '_lin'; end adiabe_ = []; if ~KIN_E; adiabe_ = '_adiabe'; end % resolution and boxsize res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)]; if (LX ~= LY) geo_ = ['_Lx_',num2str(LX),'_Ly_',num2str(LY)]; else geo_ = ['_L_',num2str(LX)]; end if (NZ > 1) %3D case res_ = [res_,'x',num2str(NZ)]; if abs(Q0) > 0 geo_ = [geo_,'_q0_',num2str(Q0)]; end if abs(EPS) > 0 geo_ = [geo_,'_e_',num2str(EPS)]; end if abs(SHEAR) > 0 geo_ = [geo_,'_s_',num2str(SHEAR)]; end end % put everything together in the param character chain u_ = '_'; % underscore variable PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_]; 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 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); 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_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.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.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,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; % system(MAKE); %% disp(['Set up ',SIMID]); disp([res_,geo_,HLdeg_]); if JOB2LOAD>=0 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 d35641f..2643945 100644 --- a/matlab/write_fort90.m +++ b/matlab/write_fort90.m @@ -1,99 +1,100 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC) % Write the input script "fort.90" with desired parameters INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.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,[' 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,[' 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,[' Npol = ', num2str(GRID.Npol),'\n']); fprintf(fid,[' SG = ', GRID.SG,'\n']); fprintf(fid,'/\n'); fprintf(fid,'&GEOMETRY\n'); fprintf(fid,[' geom = ', GEOM.geom,'\n']); fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']); fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']); fprintf(fid,[' eps = ', num2str(GEOM.eps),'\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,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']); fprintf(fid,[' write_gamma = ', OUTPUTS.write_gamma,'\n']); fprintf(fid,[' write_hf = ', OUTPUTS.write_hf,'\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,[' job2load = ', num2str(OUTPUTS.job2load),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&MODEL_PAR\n'); fprintf(fid,' ! Collisionality\n'); fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']); fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']); fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']); fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']); fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']); fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']); fprintf(fid,[' mu_z = ', num2str(MODEL.mu_z),'\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,[' K_n = ', num2str(MODEL.K_n),'\n']); fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']); fprintf(fid,[' K_E = ', num2str(MODEL.K_E),'\n']); fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']); fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']); fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&COLLISION_PAR\n'); fprintf(fid,[' collision_model = ', COLL.collision_model,'\n']); fprintf(fid,[' gyrokin_CO = ', COLL.gyrokin_CO,'\n']); fprintf(fid,[' interspecies = ', COLL.interspecies,'\n']); fprintf(fid,[' mat_file = ', COLL.mat_file,'\n']); fprintf(fid,'/\n'); fprintf(fid,'&INITIAL_CON\n'); fprintf(fid,[' INIT_OPT = ', INITIAL.INIT_OPT,'\n']); fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\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,'/\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