diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m index 47ba061..7ee54f6 100644 --- a/matlab/compute/compute_fluxtube_growth_rate.m +++ b/matlab/compute/compute_fluxtube_growth_rate.m @@ -1,73 +1,83 @@ 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; - +ikynz = (abs(DATA.PHI(:,1,1,1)) > 0); +ikynz = logical(ikynz .* (DATA.ky > 0)); phi = DATA.PHI(ikynz,ikxnz,:,:); 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,2),3); wl = log(phi_n./phi_nm1)/dt; w_ky(:,is) = squeeze(sum(sum(wl.*phi_nm1,2),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,:,:),2),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; +linear_gr.std_g = std (real(w_ky(Is,:)),[],2); +linear_gr.avg_g = mean(real(w_ky(Is,:)),2); +linear_gr.std_w = std (imag(w_ky(Is,:)),[],2); +linear_gr.avg_w = mean(imag(w_ky(Is,:)),2); + if PLOT >0 - figure + figure +if PLOT >1 + subplot(1,PLOT,1) +end + plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','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),'x','DisplayName','$\epsilon$'); hold on; + + errorbar(linear_gr.ky,linear_gr.avg_g,linear_gr.std_g,... + '-o','DisplayName','$Re(\omega_{k_y})$',... + 'LineWidth',1.5); hold on; + errorbar(linear_gr.ky,linear_gr.avg_w,linear_gr.std_w,... + '--*','DisplayName','$Im(\omega_{k_y})$',... + 'LineWidth',1.5); 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 + +if PLOT > 1 + subplot(1,PLOT,2) + plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on; + plot(DATA.Ts3D(its+1:ite),linear_gr.w_ky(Is,:)); + xlabel('t'); ylabel('$\gamma(t),\omega(t)$') +end + +if PLOT > 2 + subplot(1,PLOT,3) + semilogy(DATA.Ts3D,squeeze(abs(DATA.PHI(2,1,DATA.Nz/2,:)))); hold on; + xlabel('t'); ylabel('$|\phi_{ky}|(t)$') end end \ No newline at end of file diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m index e996e71..b9b99f3 100644 --- a/matlab/extract_fig_data.m +++ b/matlab/extract_fig_data.m @@ -1,41 +1,42 @@ fig = gcf; axObjs = fig.Children; dataObjs = axObjs.Children; X_ = []; Y_ = []; for i = 1:numel(dataObjs) X_ = [X_ dataObjs(i).XData]; Y_ = [Y_ dataObjs(i).YData]; end n0 = 1; +n1 = 26;%numel(X_); figure; mvm = @(x) movmean(x,1); shift = X_(n0); % shift = 0; % plot(X_(n0:end),Y_(n0:end)); -plot(mvm(X_(n0:end)-shift),mvm(Y_(n0:end))); hold on; +plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on; t0 = ceil(numel(X_)*0.5); t1 = numel(X_); avg= mean(Y_(t0:t1)); dev = std(Y_(t0:t1)); disp(['AVG =',sprintf('%2.2f',avg),'+-',sprintf('%2.2f',dev)]); % % n1 = n0+1; n2 = min(n1 + 50000,numel(Y_)); % avg_ = mean(Y_(n1:n2)); % std_ = std(Y_(n1:n2)); % title(['avg = ',num2str(avg_),' std = ', num2str(std_)]) % plot([X_(n1),X_(n2)],[1 1]*avg_,'--k') % ylim([0,avg_*3]); % sum(Y_(2:end)./X_(2:end).^(3/2)) %% if 0 %% m1 = mean(y1); m2 = mean(y2); s1 = std(y1); s2 = std(y2); delta = abs(m2-m1)/min(s1,s2); delta end % xlim([-3,3]); ylim([0,4.5]); \ No newline at end of file diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m index 0e0e991..3ac095c 100644 --- a/matlab/load/quick_gene_load.m +++ b/matlab/load/quick_gene_load.m @@ -1,48 +1,55 @@ %% 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.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.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'; path = '/home/ahoffman/gene/linear_CBC_results/'; -fname = 'CBC_linear.txt'; +% fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt'; +% fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt'; +% fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt'; +% fname = 'CBC_KT_5.3_20x1x32x32x12_Lv_3_Lw_12.txt'; +% fname = 'CBC_KT_5.3_32x1x48x40x16_Lv_3_Lw_12.txt'; +fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt'; +% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12_nuv_1e-3.txt'; data_ = load([path,fname]); figure -plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE'); \ No newline at end of file +plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE'); hold on; +plot(data_(:,2),data_(:,4),'--*k','DisplayName','GENE'); \ No newline at end of file diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m index 14ebd49..7288fe1 100644 --- a/matlab/plot/spectrum_1D.m +++ b/matlab/plot/spectrum_1D.m @@ -1,164 +1,167 @@ function [ FIGURE ] = spectrum_1D( data, options ) options.PLAN = 'kxky'; options.COMP = 'avg'; options.INTERP = 0; options.POLARPLOT = 0; options.AXISEQUAL = 1; switch options.COMPT case 'avg' [~,it0] = min(abs(data.Ts3D-options.TIME(1))); [~,it1] = min(abs(data.Ts3D-options.TIME(end))); options.TIME = data.Ts3D(it0:it1); end toplot = process_field(data,options); t = data.Ts3D; frames = toplot.FRAMES; colors = jet(numel(frames)); switch options.NAME case '\Gamma_x' FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS]; fieldname = 'transport'; case '\phi' FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS]; fieldname = 'ES pot.'; case 'n_i' FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS]; yname = 'n_i'; fieldname = 'ion dens.'; case 'n_e' FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS]; fieldname = 'elec. dens.'; case 'N_i^{00}' FIGURE.fig = figure; FIGURE.FIGNAME = ['Ni00_spectrum_',data.PARAMS]; fieldname = 'ion gdens.'; end PLOT2D = 0; switch options.COMPXY case 'avg' compx = @(x) mean(x,2); compy = @(x) mean(x,1); ynamecx = ['$\sum_{k_x}',options.NAME,'$']; ynamecy = ['$\sum_{k_y}',options.NAME,'$']; case 'sum' compx = @(x) sum(x,2); compy = @(x) sum(x,1); ynamecx = ['$\sum_{k_x}',options.NAME,'$']; ynamecy = ['$\sum_{k_y}',options.NAME,'$']; case 'max' compx = @(x) max(x,2); compy = @(x) max(x,1); ynamecx = ['$\max_{k_x}',options.NAME,'$']; ynamecy = ['$\max_{k_y}',options.NAME,'$']; case 'zero' compx = @(x) x(:,data.Nx/2); compy = @(x) x(1,:); ynamecx = ['$',options.NAME,'(k_x=',num2str(data.kx(1)),')$']; ynamecy = ['$',options.NAME,'(k_y=',num2str(data.ky(1)),')$']; otherwise compx = @(x) x(:,:); compy = @(x) x(:,:); PLOT2D= 1; end if ~PLOT2D set(gcf, 'Position', [20 50 1200 500]) subplot(1,2,1) k = data.ky; xname = '$k_y$'; nmax = ceil(data.Nkx*2/3); shiftx = @(x) x;%(1:nmax); shifty = @(x) x;%(1:nmax); switch options.COMPT case 'avg' it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); Gk = compx(abs(mean(toplot.FIELD(:,:,:),3))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = shiftx(k); if options.OK Y = shifty(Gk)./X; else Y = shifty(Gk); end plot(X,Y,'DisplayName','t-averaged') otherwise for it = 1:numel(frames) Gk = compx(abs(toplot.FIELD(:,:,it))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = shiftx(k); if options.OK Y = shifty(Gk)./X; else Y = shifty(Gk); end plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],... 'Color',colors(it,:)); hold on; end end grid on title(['GM $k_x$ ',fieldname]); % legend('show','Location','eastoutside') xlabel(xname); ylabel(ynamecx) subplot(1,2,2) k = data.kx; xname = '$k_x$'; nmax = floor(data.Nky/2*2/3); switch options.COMPT case 'avg' % it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); Gk = compy(abs(mean(toplot.FIELD(:,:,:),3)))'; Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = k(k>0); X = X(1:end-1); Yp= Gk(k>0); Ym= Gk(k<0); Y = Yp(1:end-1) + Ym(end:-1:1); Y = Y(end:-1:1); plot(X,Y,'DisplayName','t-averaged') otherwise for it = 1:numel(toplot.FRAMES) Gk = compy(abs(toplot.FIELD(:,:,it))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = k(k>0); X = X(1:end-1); Yp= Gk(k>0); Ym= Gk(k<0); Y = Yp(1:end-1) + Ym(end:-1:1); Y = Y(end:-1:1); plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],... 'Color',colors(it,:)); hold on; end end grid on title(['GM $k_y$ ',fieldname]); % legend('show','Location','eastoutside'); xlabel(xname); ylabel(ynamecy) else % it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); Gk = mean(abs(toplot.FIELD(:,:,:)),3); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end pclr = pcolor(toplot.X,toplot.Y,Gk); set(pclr, 'edgecolor','none'); + xlabel('$k_x$'); + ylabel('$k_y$'); + title(['GM $|$',fieldname,'$(k_x,k_y)|$ t-averaged']); end diff --git a/matlab/setup.m b/matlab/setup.m index 6228cac..db703d0 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,177 +1,177 @@ %% _______________________________________________________________________ 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.beta = BETA; MODEL.mu_x = MU_X; MODEL.mu_y = MU_Y; MODEL.N_HD = N_HD; MODEL.mu_z = MU_Z; 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.K_E = 0; % 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_ = ['_nu_%1.1e_',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/wk/CBC_kT_PJ_scan_quick_run.m b/wk/CBC_kT_PJ_scan_quick_run.m new file mode 100644 index 0000000..4245d7d --- /dev/null +++ b/wk/CBC_kT_PJ_scan_quick_run.m @@ -0,0 +1,47 @@ +KT_a = [3:0.5:7]; +g_max= KT_a*0; +g_avg= KT_a*0; +g_std= KT_a*0; +k_max= KT_a*0; + +NU = 0.05; +DT = 2e-3; +TMAX = 50; +ky_ = 0.3; +SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation +% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation +RUN = 1; +figure + +for P = [20] + +i=1; +for K_T = KT_a + + quick_run + + g_max(i) = gmax; + k_max(i) = kmax; + + g_avg(i) = lg.avg_g; + g_std(i) = lg.std_g; + + i = i + 1; +end +%% + +% plot(KT_a,max(g_max,0)); +y_ = g_avg; +e_ = g_std; + +y_ = y_.*(y_-e_>0); +e_ = e_ .* (y_>0); +errorbar(KT_a,y_,e_,... + 'LineWidth',1.2,... + 'DisplayName',['(',num2str(P),',',num2str(P/2),')']); +hold on; +title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']); +legend('show'); xlabel('$K_T$'); ylabel('$\gamma$'); +drawnow +end + diff --git a/wk/CBC_nu_CO_scan_quick_run.m b/wk/CBC_nu_CO_scan_quick_run.m new file mode 100644 index 0000000..6b9bc4f --- /dev/null +++ b/wk/CBC_nu_CO_scan_quick_run.m @@ -0,0 +1,48 @@ +% NU_a = [0.05 0.15 0.25 0.35 0.45]; +NU_a = [0:0.05:0.5]; +g_max= NU_a*0; +g_avg= NU_a*0; +g_std= NU_a*0; +k_max= NU_a*0; +CO = 'LR'; + +K_T = 5.3; +DT = 2e-3; +TMAX = 30; +ky_ = 0.3; +SIMID = 'linear_CBC_nu_scan_ky=0.3_CLOS_0_LRGK'; % Name of the simulation +RUN = 1; +figure + +for P = [2 4 6 10 12 20] + +i=1; +for NU = NU_a + + quick_run + + g_max(i) = gmax; + k_max(i) = kmax; + + g_avg(i) = lg.avg_g; + g_std(i) = lg.std_g; + + i = i + 1; +end +%% + +% plot(KT_a,max(g_max,0)); +y_ = g_avg; +e_ = g_std; + +y_ = y_.*(y_-e_>0); +e_ = e_ .* (y_>0); +errorbar(NU_a,y_,e_,... + 'LineWidth',1.2,... + 'DisplayName',['(',num2str(P),',',num2str(P/2),')']); +hold on; +title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']); +legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$'); +drawnow +end + diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m new file mode 100644 index 0000000..96c8b4b --- /dev/null +++ b/wk/Dimits_2000_fig3.m @@ -0,0 +1,28 @@ +%192x96x16x3x2 simulation +kT_Xi_GM_32 = ... + [7.0 3.6;... + 6.0 2.6;... + 5.0 1.1;... + ]; +%128x64x16x5x3 simulation +kT_Xi_GM_53 = ... + [7.0 4.6;... + 5.3 1.9;... + 5.0 1.5;... + ]; +%128x64x16x24x12 simulation +kT_Xi_GN = ... + [7.0 5.0;... + 5.3 2.4;... + 4.5 nan;... + ]; + +figure; +plot(kT_Xi_GM_32(:,1), kT_Xi_GM_32(:,2),'o-','DisplayName','192x96x16x3x2'); hold on +plot(kT_Xi_GM_53(:,1), kT_Xi_GM_53(:,2),'o-','DisplayName','128x64x16x5x3'); hold on +plot(kT_Xi_GN(:,1), kT_Xi_GN(:,2),'o--k','DisplayName','GENE 128x64x16x24x12'); +xlabel('$\kappa_T$'); ylabel('$\chi_i$'); +xlim([0 20]); +ylim([0 12]); +title('Dimits et al. 2000, Fig. 3'); +legend('show'); grid on; \ No newline at end of file diff --git a/wk/ITG_peak_KT_53_k_035_scan.m b/wk/ITG_peak_KT_53_k_035_scan.m new file mode 100644 index 0000000..9671fc0 --- /dev/null +++ b/wk/ITG_peak_KT_53_k_035_scan.m @@ -0,0 +1,6 @@ +P = [ 2 4 8 10 16]; + +k = 0.35; KN = 2.22; KT = 5.3; + +g = [0.22 0.29 0.10 0.04 0.11]; +w = [0.00 0.00 0.00 0.00 0.00]; \ No newline at end of file diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m index 7721a36..1ba1f5a 100644 --- a/wk/analysis_HeLaZ.m +++ b/wk/analysis_HeLaZ.m @@ -1,207 +1,207 @@ addpath(genpath([helazdir,'matlab'])) % ... add addpath(genpath([helazdir,'matlab/plot'])) % ... add addpath(genpath([helazdir,'matlab/compute'])) % ... add addpath(genpath([helazdir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [helazdir,'results/',outfile,'/']; MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax -JOBNUMMIN = 00; JOBNUMMAX = 20; +JOBNUMMIN = 00; JOBNUMMAX = 10; data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) -options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%(1/data.Nx/data.Ny)^2; +options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2; options.TAVG_1 = data.Ts3D(end); % Averaging times duration options.NMVA = 1; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 1; fig = plot_radial_transport_and_spacetime(data,options); save_figure(data,fig,'.png') end if 0 %% statistical transport averaging options.T = [200 400]; fig = statistical_transport_averaging(data,options); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xz'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; % options.TIME = data.Ts5D(end-30:end); -options.TIME = data.Ts3D(1:2:end); -% options.TIME = [750:1:1000]; +% options.TIME = data.Ts3D(1:10:end); +options.TIME = [200:2:400]; data.EPS = 0.1; data.a = data.EPS * 2000; create_film(data,options,'.gif') end if 0 %% 2D snapshots % Options options.INTERP = 0; options.POLARPLOT = 0; -options.AXISEQUAL = 1; +options.AXISEQUAL = 0; % options.NAME = '\phi'; % options.NAME = 'n_e'; options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; -options.PLAN = 'kxz'; +options.PLAN = 'kxky'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; -options.TIME = [120 150 155]; +options.TIME = [20 60 200 400 500]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 1; options.NAME = '\phi'; options.PLANES = [1:1:16]; options.TIME = [15]; options.PLT_MTOPO = 0; data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) -options.SPAR = linspace(-3,3,32)+(6/127/2); -options.XPERP = linspace( 0,6,32); -% options.SPAR = gene_data.vp'; -% options.XPERP = gene_data.mu'; -options.iz = 'avg'; -options.T = [600]; +% options.SPAR = linspace(-3,3,32)+(6/127/2); +% options.XPERP = linspace( 0,6,32); +options.SPAR = gene_data.vp'; +options.XPERP = gene_data.mu'; +options.iz = 1; +options.T = [500 1000]; options.PLT_FCT = 'contour'; options.ONED = 0; options.non_adiab = 0; options.SPECIE = 'i'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; -options.NORMALIZED = 1; +options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [1000]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum -options.TIME = [300 600]; +options.TIME = [100 500]; options.NORM =1; % options.NAME = '\phi'; -% options.NAME = 'N_i^{00}'; -options.NAME ='\Gamma_x'; +options.NAME = 'N_i^{00}'; +% options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = '2D'; % avg/sum/max/zero/ 2D plot otherwise options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 'avg'; options.COMPY = 'avg'; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:160]; options.NMA = 1; -options.NMODES = 15; +options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos fig = plot_metric(data); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 15]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig,'.png') end diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m index 4b20e0b..57ba009 100644 --- a/wk/analysis_gene.m +++ b/wk/analysis_gene.m @@ -1,142 +1,145 @@ % folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/'; % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/rm_corrections_HF/'; % folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/'; % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/'; % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/'; % folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/'; % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/'; % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/'; % folder = '/misc/gene_results/LD_zpinch_1.6/'; % folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/'; % folder = '/misc/gene_results/ZP_kn_2.5_large_box/'; folder = '/misc/gene_results/CBC/128x64x16x24x12/'; +% folder = '/misc/gene_results/CBC/196x96x20x32x16_00/'; +% folder = '/misc/gene_results/CBC/128x64x16x6x4/'; +% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/'; gene_data = load_gene_data(folder); gene_data = invert_kxky_to_kykx_gene_results(gene_data); if 1 %% Space time diagramm (fig 11 Ivanov 2020) options.TAVG_0 = 0.6*gene_data.Ts3D(end); options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration options.NMVA = 1; % Moving average for time traces options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x) options.INTERP = 1; gene_data.FIGDIR = folder; fig = plot_radial_transport_and_spacetime(gene_data,options); save_figure(gene_data,fig,'.png') end if 0 %% statistical transport averaging options.T = [100 500]; fig = statistical_transport_averaging(gene_data,options); end if 0 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; % options.NAME = 'Q_x'; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'xy'; % options.NAME ='f_e'; % options.PLAN = 'sx'; options.COMP = 1; options.TIME = [15 50 100 200]; gene_data.a = data.EPS * 2000; fig = photomaton(gene_data,options); save_figure(gene_data,fig,'.png') end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = gene_data.Ts3D; gene_data.a = data.EPS * 2000; create_film(gene_data,options,'.gif') end if 0 %% Geometry names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',... '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',... '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'}; figure; subplot(311) for i = 1:6 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry'); subplot(312) for i = 7:10 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); subplot(313) for i = 11:16 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); end if 0 %% Show f_i(vpar,mu) -options.times = 200:300; +options.times = 20:80; options.specie = 'i'; options.PLT_FCT = 'contour'; options.folder = folder; options.iz = 1; options.FIELD = ''; options.ONED = 0; % options.FIELD = 'Q_es'; plot_fa_gene(options); end if 0 %% Time averaged spectrum options.TIME = 300:600; options.NORM =1; % options.NAME = '\phi'; % options.NAME = 'n_i'; options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'zero'; options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(gene_data,options); % save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = 100:700; options.NMA = 1; options.NMODES = 15; options.iz = 'avg'; fig = mode_growth_meter(gene_data,options); save_figure(gene_data,fig) end diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index 9b6139f..45ec1e0 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,45 +1,52 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" helazdir = '/home/ahoffman/HeLaZ/'; % Directory of the simulation (from results) % if 1% Local results % outfile ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits % outfile ='shearless_cyclone/128x64x16x5x3_Dim_90'; % outfile ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; % outfile ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; % outfile ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS % outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; % outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; % outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; % outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; %% shearless CBC % outfile ='shearless_cyclone/64x32x16x5x3_CBC_080'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_120'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_080'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_100'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_120'; % outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% ZPINCH % outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3'; % outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3'; % outfile ='Zpinch_rerun/Kn_2.5_312x196x5x3_Lx_400_Ly_200'; % outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3'; -outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box'; +% outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box'; % outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120'; % outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4'; % outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6'; % outfile ='Zpinch_rerun/Kn_1.6_256x128x21x11'; %% CBC % outfile = 'CBC/64x32x16x5x3'; % outfile = 'CBC/64x128x16x5x3'; % outfile = 'CBC/128x64x16x5x3'; -outfile = 'CBC/64x64x16x3x2'; +% outfile = 'CBC/192x96x16x3x2'; +% outfile = 'CBC/128x64x16x5x3'; + +outfile = 'CBC/kT_scan_128x64x16x5x3'; +% outfile = 'CBC/kT_scan_192x96x16x3x2'; + +%% Linear CBC +% outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; run analysis_HeLaZ diff --git a/wk/quick_run.m b/wk/quick_run.m index e0bcee2..8a7931f 100644 --- a/wk/quick_run.m +++ b/wk/quick_run.m @@ -1,183 +1,187 @@ %% QUICK RUN SCRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % -RUN = 1; % To run or just to load +% RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options HELAZDIR = '/home/ahoffman/HeLaZ/'; EXECNAME = 'helaz3'; -% EXECNAME = 'helaz3_shear'; %verified version +% EXECNAME = 'helaz3_dbg'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.01; % Collision frequency +% NU = 0.0; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_N = 2.22;%2.0; % Density gradient drive -K_T = 6.92;%0.25*K_N; % Temperature ''' -K_E = 0.0; % Electrostat ''' +K_N = 2.22;%2.0; % ion Density gradient drive +K_Ne = K_N; % ele Density gradient drive +% K_T = 4.0;%0.25*K_N; % Temperature ''' +K_Te = K_T; % Temperature ''' % SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0e-1; % electron plasma beta %% GRID PARAMETERS -PMAXE = 4; % Hermite basis size of electrons -JMAXE = 2; % Laguerre " -PMAXI = 4; % " ions -JMAXI = 2; % " +% P = 12; +J = P/2; +PMAXE = P; % Hermite basis size of electrons +JMAXE = J; % Laguerre " +PMAXI = P; % " ions +JMAXI = J; % " NX = 12; % real space x-gridpoints -NY = 8; % '' y-gridpoints +NY = 2; % '' y-gridpoints LX = 2*pi/0.1; % Size of the squared frequency domain -LY = 2*pi/0.1; % Size of the squared frequency domain +% LY = 2*pi/0.3; % Size of the squared frequency domain +LY = 2*pi/ky_; NZ = 16; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY % GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps GEOMETRY= 's-alpha'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear (Not implemented yet) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 20; % Maximal time unit -DT = 5e-3; % Time step +% TMAX = 100; % Maximal time unit +% DT = 2e-3; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 1; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS -SIMID = 'linear_CBC'; % Name of the simulation +% SIMID = 'linear_CBC'; % Name of the simulation +% SIMID = 'scan'; % Name of the simulation LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) -CO = 'DG'; -GKCO = 0; % gyrokinetic operator +% CO = 'SG'; +GKCO = 1; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; -CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s +CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio kmax = NX*pi/LX;% Highest fourier mode MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN % system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) - system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) + system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [HELAZDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) trange = [0.5 1]*data.Ts3D(end); -nplots = 1; +nplots = 0; lg = compute_fluxtube_growth_rate(data,trange,nplots); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 0 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.1 0.2 0.3]; options.normalized = 1; options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end - if 0 %% Mode evolution -options.NORMALIZED = 1; +options.NORMALIZED = 0; options.K2PLOT = 1; -options.TIME = [0 1]*data.Ts3D(end); +options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; -options.iz = 9; +options.iz = 'avg'; fig = mode_growth_meter(data,options); -save_figure(gbms_dat,fig) +save_figure(data,fig,'.png') end if 0 %% RH TEST ikx = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx))) end \ No newline at end of file diff --git a/wk/scan_quick_run.m b/wk/scan_quick_run.m new file mode 100644 index 0000000..4245d7d --- /dev/null +++ b/wk/scan_quick_run.m @@ -0,0 +1,47 @@ +KT_a = [3:0.5:7]; +g_max= KT_a*0; +g_avg= KT_a*0; +g_std= KT_a*0; +k_max= KT_a*0; + +NU = 0.05; +DT = 2e-3; +TMAX = 50; +ky_ = 0.3; +SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation +% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation +RUN = 1; +figure + +for P = [20] + +i=1; +for K_T = KT_a + + quick_run + + g_max(i) = gmax; + k_max(i) = kmax; + + g_avg(i) = lg.avg_g; + g_std(i) = lg.std_g; + + i = i + 1; +end +%% + +% plot(KT_a,max(g_max,0)); +y_ = g_avg; +e_ = g_std; + +y_ = y_.*(y_-e_>0); +e_ = e_ .* (y_>0); +errorbar(KT_a,y_,e_,... + 'LineWidth',1.2,... + 'DisplayName',['(',num2str(P),',',num2str(P/2),')']); +hold on; +title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']); +legend('show'); xlabel('$K_T$'); ylabel('$\gamma$'); +drawnow +end +