diff --git a/matlab/compile_results.m b/matlab/compile_results.m index 3e4bb91..a2689bc 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,89 +1,103 @@ CONTINUE = 1; JOBNUM = 0; JOBFOUND = 0; +TJOB_SE = []; % Start and end times of jobs +NU_EVOL = []; % evolution of parameter nu between jobs +ETAB_EVOL= []; % +L_EVOL = []; % + +% FIELDS Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; GGAMMA_ = []; PGAMMA_ = []; PHI_ = []; Ts0D_ = []; Ts2D_ = []; Ts5D_ = []; Sipj_ = []; Sepj_ = []; Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old; while(CONTINUE) filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); if exist(filename, 'file') == 2 % Load results of simulation #JOBNUM load_results % Check polynomials degrees Pe_new= numel(Pe); Je_new= numel(Je); Pi_new= numel(Pi); Ji_new= numel(Ji); % If 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 end if W_GAMMA GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); Ts0D_ = cat(1,Ts0D_,Ts0D); end if W_PHI || W_NA00 Ts2D_ = cat(1,Ts2D_,Ts2D); end if W_PHI PHI_ = cat(3,PHI_,PHI); end if W_NA00 Ni00_ = cat(3,Ni00_,Ni00); Ne00_ = cat(3,Ne00_,Ne00); end if W_NAPJ || W_SAPJ Ts5D_ = cat(1,Ts5D_,Ts5D); end if W_NAPJ Nipj_ = cat(5,Nipj_,Nipj); Nepj_ = cat(5,Nepj_,Nepj); end if W_SAPJ Sipj_ = cat(5,Sipj_,Sipj); Sepj_ = cat(5,Sepj_,Sepj); end + % Evolution of simulation parameters + load_params + TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)]; + NU_EVOL = [NU_EVOL NU NU]; + ETAB_EVOL = [ETAB_EVOL ETAB ETAB]; + L_EVOL = [L_EVOL L L]; + JOBFOUND = JOBFOUND + 1; LASTJOB = JOBNUM; elseif (JOBNUM > 20) CONTINUE = 0; disp(['found ',num2str(JOBFOUND),' results']); end JOBNUM = JOBNUM + 1; Pe_old = Pe_new; Je_old = Je_new; Pi_old = Pi_new; Ji_old = Ji_new; + end GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_; Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_; Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_; clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_ Sipj = Sipj_; Sepj = Sepj_; clear Sipj_ Sepj_ JOBNUM = LASTJOB filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); \ No newline at end of file diff --git a/matlab/create_gif_1D_phi.m b/matlab/create_gif_1D_phi.m new file mode 100644 index 0000000..f68761c --- /dev/null +++ b/matlab/create_gif_1D_phi.m @@ -0,0 +1,62 @@ +title1 = GIFNAME; +FIGDIR = BASIC.RESDIR; +if ~exist(FIGDIR, 'dir') + mkdir(FIGDIR) +end + +GIFNAME = [FIGDIR, GIFNAME,'.gif']; + +% Set colormap boundaries +hmax = max(max(max(FIELD))); +hmin = min(min(min(FIELD))); +fieldabsmin = min(min(abs(FIELD(:,:)))); +fieldabsmax = max(max(abs(FIELD(:,:)))); +flag = 0; +if hmax == hmin + disp('Warning : h = hmin = hmax = const') +else +% Setup figure frame +fig = figure; +subplot(211) + plot(X,FIELD(:,1)); % to set up + axis tight manual % this ensures that getframe() returns a consistent size + in = 1; + nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); + for n = FRAMES % loop over selected frames + scale = max(FIELD(:,n))*SCALING + (1-SCALING); + + subplot(211) + plot(X,FIELD(:,n)/scale,linestyle); + if (YMIN ~= YMAX && XMIN ~= XMAX) + ylim([YMIN,YMAX]); xlim([XMIN,XMAX]); + end + title(['$t \approx$', sprintf('%.3d',ceil(T(n))), ', scaling = ',sprintf('%.1e',scale)]); + xlabel(XNAME); ylabel(FIELDNAME); + drawnow + % Capture the plot as an image + frame = getframe(fig); + im = frame2im(frame); + [imind,cm] = rgb2ind(im,64); + % Write to the GIF File + if in == 1 + imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf); + else + imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY); + end + % terminal info + while nbytes > 0 + fprintf('\b') + nbytes = nbytes - 1; + end + nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:))); + in = in + 1; + + subplot(212) + ylim([fieldabsmin,fieldabsmax]); xlim([T(FRAMES(1)),T(FRAMES(end))]); hold on; + plot(T(n),abs(max(FIELD(:,n))),'.k'); xlabel('$\rho_s/c_s$') + + end + disp(' ') + disp(['Gif saved @ : ',GIFNAME]) +end + diff --git a/matlab/profiler.m b/matlab/profiler.m index 37d363b..57b8470 100644 --- a/matlab/profiler.m +++ b/matlab/profiler.m @@ -1,61 +1,95 @@ %% load profiling % filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00); CPUTIME = double(h5readatt(filename,'/data/input','cpu_time')); DT_SIM = h5readatt(filename,'/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'); missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ... -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_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]); 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)); +diag_Ta = mean(diff(diag_Tc)); + +NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM; + %% Plots -TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc); diff(ghost_Tc);... - diff(Sapj_Tc); diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)]; -TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/8,8]); +if 1 + %% Area plot fig = figure; p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none'); -legend('Compute RHS','Adv. fields','Poisson', 'comm', 'Sapj','Check+Sym','Diag','Missing') +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') 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 -%% Plots -% fig = figure; -% -% p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none'); -% legend('Compute RHS','Adv. fields','Poisson','Sapj','Check+Sym','Diag','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 \ No newline at end of file +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; +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 +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(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') +xlabel('Step Comp. Time [s]'); ylabel('') +FIGNAME = 'profiler'; +save_figure +end \ No newline at end of file diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m index 0425cf4..98c1e35 100644 --- a/matlab/write_sbash_daint.m +++ b/matlab/write_sbash_daint.m @@ -1,56 +1,57 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $SCRATCH/HeLaZ/wk\n',... ... 'cd $SCRATCH/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... 'cp $HOME/HeLaZ/wk/fort.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... 'jid=$(sbatch batch_script.sh)\n',... 'echo $jid\n',... 'echo to check output log :\n',... 'echo tail -f $SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --job-name="',CLUSTER.JNAME,'"\n',... '#SBATCH --time=', CLUSTER.TIME,'\n',... '#SBATCH --nodes=', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',... '#SBATCH --ntasks-per-core=', CLUSTER.NTPC,'\n',... '#SBATCH --mem=', CLUSTER.MEM,'\n',... '#SBATCH --error=err.txt\n',... '#SBATCH --output=out.txt\n',... '#SBATCH --account="s882"\n',... '#SBATCH --constraint=mc\n',... '#SBATCH --hint=nomultithread\n',... '#SBATCH --partition=',CLUSTER.PART,'\n',... 'export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n',... ...% '#SBATCH --job-name=',PARAMS,'\n\n',... 'module purge\n',... 'module load PrgEnv-intel\n',... 'module load cray-hdf5-parallel\n',... 'module load cray-mpich\n',... 'module load craype-x86-skylake\n',... 'module load cray-fftw\n',... -'srun ./../../../bin/helaz']); +'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]); +%'srun ./../../../bin/helaz']); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk'); \ No newline at end of file diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m index dd7db86..b4fdf01 100644 --- a/matlab/write_sbash_marconi.m +++ b/matlab/write_sbash_marconi.m @@ -1,47 +1,44 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $CINECA_SCRATCH/HeLaZ/wk\n',... ... 'cd $CINECA_SCRATCH/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... 'cp $HOME/HeLaZ/wk/fort.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... 'sbatch batch_script.sh\n',... 'echo tail -f $CINECA_SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --job-name=',CLUSTER.JNAME,'\n',... '#SBATCH --time=', CLUSTER.TIME,'\n',... '#SBATCH --nodes=', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',... '#SBATCH --mem=', CLUSTER.MEM,'\n',... '#SBATCH --error=err.txt\n',... '#SBATCH --output=out.txt\n',... '#SBATCH --account=FUA35_TSVVT421\n',... '#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',... -'module load intel\n',... -'module load intelmpi\n',... -'module load autoload hdf5/1.10.4--intelmpi--2018--binary\n',... -'module load fftw\n',... +'module load autoload hdf5 fftw\n',... 'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk'); \ No newline at end of file diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index 50769e0..4eca5ea 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,454 +1,466 @@ %% Load results outfile =''; -if 1 +if 0 %% Load from Marconi -outfile =''; -% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_20_J_3_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; -% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.8_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.8_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; -% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.7_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.7_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; -outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; -outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/out.txt'; +outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_eta_0.6_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; +% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_80_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_4e-03/out.txt'; +% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e+00/200x100_L_120_P_12_J_6_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_2e-02/out.txt'; +% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; +% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/out.txt'; +% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; BASIC.RESDIR = load_marconi(outfile); end if 0 %% Load from Daint - outfile =''; + outfile ='/scratch/snx3000/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt'; BASIC.RESDIR = load_daint(outfile); end %% -% JOBNUM = 0; load_results; % JOBNUM = 1; load_results; compile_results -load_params %% Retrieving max polynomial degree and sampling info Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe); Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi); Ns5D = numel(Ts5D); Ns2D = numel(Ts2D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts2D = Ts2D'; %% Build grids Nkr = numel(kr); Nkz = numel(kz); [KZ,KR] = meshgrid(kz,kr); Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz); dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1); KPERP2 = KZ.^2+KR.^2; [~,ikr0] = min(abs(kr)); [~,ikz0] = min(abs(kz)); Lk = max(Lkr,Lkz); dr = 2*pi/Lk; dz = 2*pi/Lk; Nr = max(Nkr,Nkz); Nz = Nr; r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r); z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z); [ZZ,RR] = meshgrid(z,r); %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nr,Nz,Ns2D); % Gyrocenter density ni00 = zeros(Nr,Nz,Ns2D); np_i = zeros(Nr,Nz,Ns5D); % Ion particle density si00 = zeros(Nr,Nz,Ns5D); phi = zeros(Nr,Nz,Ns2D); drphi = zeros(Nr,Nz,Ns2D); dr2phi = zeros(Nr,Nz,Ns2D); dzphi = zeros(Nr,Nz,Ns2D); for it = 1:numel(Ts2D) NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz))); ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz))); phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz))); drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz))); dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz))); dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz))); end % Building a version of phi only 5D sampling times PHI_Ts5D = zeros(Nkr,Nkz,Ns5D); err = 0; for it = 1:numel(Ts5D) % Loop over 5D arrays [shift, it2D] = min(abs(Ts2D-Ts5D(it))); if shift > err; err = shift; end; PHI_Ts5D(:,:,it) = PHI(:,:,it2D); end if err > 0; disp('WARNING Ts2D and Ts5D are shifted'); end; Np_i = zeros(Nkr,Nkz,Ns5D); % Ion particle density in Fourier space for it = 1:numel(Ts5D) [~, it2D] = min(abs(Ts2D-Ts5D(it))); Np_i(:,:,it) = 0; for ij = 1:Nji Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1)); Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it)); end np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nr,Nz))); end % Post processing disp('- post processing') E_pot = zeros(1,Ns2D); % Potential energy n^2 E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2 ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi| % gyrocenter and particle flux from real space GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_re = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_ze = zeros(1,Ns2D); % Gyrocenter flux Gamma = PFlux_ri = zeros(1,Ns5D); % Particle flux % gyrocenter and particle flux from fourier coefficients GFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Ni00.*conj(PHI),1),2)))*(2*pi/Nr/Nz)^2; PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz)^2; % Hermite energy spectrum epsilon_e_pj = zeros(Npe,Nje,Ns5D); epsilon_i_pj = zeros(Npi,Nji,Ns5D); -phi_max = zeros(1,Ns2D); % Time evol. of the norm of phi +phi_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi +phi_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi +phi_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi +phi_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi Ne_norm = zeros(Npe,Nje,Ns5D); % Time evol. of the norm of Napj Ni_norm = zeros(Npi,Nji,Ns5D); % . Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2; for it = 1:numel(Ts2D) % Loop over 2D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); - phi_max(it) = max(max(squeeze(phi(:,:,it)))); + phi_maxr_maxz(it) = max( max(squeeze(phi(:,:,it)))); + phi_avgr_maxz(it) = max(mean(squeeze(phi(:,:,it)))); + phi_maxr_avgz(it) = mean( max(squeeze(phi(:,:,it)))); + phi_avgr_avgz(it) = mean(mean(squeeze(phi(:,:,it)))); ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz)))); GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; end for it = 1:numel(Ts5D) % Loop over 5D arrays [~, it2D] = min(abs(Ts2D-Ts5D(it))); Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz; epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4); epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4); % Particle flux PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lr/Lz; end %% Compute growth rate disp('- growth rate') % Find max value of transport (end of linear mode) -[~,itmax] = max(GFlux_ri); -tstart = 0.1 * Ts2D(itmax); tend = 0.9 * Ts2D(itmax); +[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nr/Nz)^2); +[~,itmax] = min(abs(Ts2D-tmax)); +tstart = 0.1 * Ts2D(itmax); tend = 0.5 * Ts2D(itmax); g_ = zeros(Nkr,Nkz); for ikr = 1:Nkr for ikz = 1:Nkz g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend); end end [gmax,ikzmax] = max(g_(1,:)); kzmax = abs(kz(ikzmax)); Bohm_transport = ETAB/ETAN*gmax/kzmax^2; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Time evolutions and growth rate fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS]; set(gcf, 'Position', [100, 100, 900, 800]) +subplot(111); + suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),... + ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$']); subplot(421); for ip = 1:Npe for ij = 1:Nje plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$']; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,... 'Color',clr,'LineStyle',lstyle{1}); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); subplot(423) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$']; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,... 'Color',clr,'LineStyle',lstyle{1}); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$') subplot(222) plot(Ts0D,GGAMMA_RI*(2*pi/Nr/Nz)^2); hold on; % plot(Ts2D,GFLUX_RI) plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2); % plot(Ts5D,PFLUX_RI,'--'); legend(['Gyro. flux';'Part. flux']); grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$') +% ylim([0,2.0]); subplot(223) - plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on; - grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show'); + plot(kz,g_(1,:),'-','DisplayName','Linear growth rate'); hold on; + plot([max(kz)*2/3,max(kz)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA'); + grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show'); +% ylim([0,max(g_(1,:))]); xlim([0,max(kz)]); subplot(224) - plotname = '$\max_{r,z}(\phi)(t)$'; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); - plot(Ts2D,phi_max,'DisplayName',plotname); hold on; - grid on; xlabel('$t c_s/R$'); ylabel('$\max_{r,z}(\phi)$'); %legend('show'); -suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]); + plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on; + plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on; + plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on; + plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on; + grid on; xlabel('$t c_s/R$'); ylabel('$T_e/e$'); legend('show'); save_figure end if 0 %% Particle fluxes SCALING = Nkr*dkr * Nkz*dkz; fig = figure; FIGNAME = ['gamma',sprintf('_%.2d',JOBNUM),'_',PARAMS]; set(gcf, 'Position', [100, 100, 800, 300]) semilogy(Ts2D,GFLUX_RI, 'color', line_colors(2,:)); hold on plot(Ts5D,PFLUX_RI,'.', 'color', line_colors(2,:)); hold on plot(Ts2D,SCALING*GFlux_ri, 'color', line_colors(1,:)); hold on plot(Ts5D,SCALING*PFlux_ri,'.', 'color', line_colors(1,:)); hold on xlabel('$tc_{s0}/\rho_s$'); ylabel('$\Gamma_r$'); grid on title(['$\eta=',num2str(ETAB),'\quad',... '\nu_{',CONAME,'}=',num2str(NU),'$', ' $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) legend('Gyro Flux','Particle flux', 'iFFT GFlux', 'iFFT PFlux')%'$\eta\gamma_{max}/k_{max}^2$') save_figure end if 0 %% Space time diagramm (fig 11 Ivanov 2020) fig = figure; FIGNAME = ['space_time_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) subplot(311) semilogy(Ts2D,GFLUX_RI,'-'); hold on plot(Ts5D,PFLUX_RI,'.'); hold on % plot(Ts2D,Bohm_transport*ones(size(Ts2D)),'--'); hold on ylabel('$\Gamma_r$'); grid on title(['$\eta=',num2str(ETAB),'\quad',... '\nu_{',CONAME,'}=',num2str(NU),'$']) legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'],'Particle flux')%'$\eta\gamma_{max}/k_{max}^2$') % set(gca,'xticklabel',[]) subplot(312) yyaxis left plot(Ts2D,squeeze(max(max((phi))))) ylabel('$\max \phi$') yyaxis right plot(Ts2D,squeeze(mean(max(dr2phi)))) ylabel('$s\sim\langle\partial_r^2\phi\rangle_z$'); grid on % set(gca,'xticklabel',[]) subplot(313) [TY,TX] = meshgrid(r,Ts2D); pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); %colorbar; xlabel('$t c_s/R$'), ylabel('$r/\rho_s$') legend('$\langle\partial_r \phi\rangle_z$') save_figure end if 0 %% Photomaton : real space % FIELD = ni00; FNAME = 'ni'; XX = RR; YY = ZZ; FIELD = phi; FNAME = 'phi'; XX = RR; YY = ZZ; % FIELD = fftshift(abs(Ni00),2); FNAME = 'Fni'; XX = fftshift(KR,2); YY = fftshift(KZ,2); % FIELD = fftshift(abs(PHI),2); FNAME = 'Fphi'; XX = fftshift(KR,2); YY = fftshift(KZ,2); tf = 100; [~,it1] = min(abs(Ts2D-tf)); tf = 118; [~,it2] = min(abs(Ts2D-tf)); tf = 140; [~,it3] = min(abs(Ts2D-tf)); tf = 300; [~,it4] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 400]) plt = @(x) x;%./max(max(x)); subplot(141) DATA = plt(FIELD(:,:,it1)); pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap gray xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it1))); subplot(142) DATA = plt(FIELD(:,:,it2)); pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap gray xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); subplot(143) DATA = plt(FIELD(:,:,it3)); pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap gray xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it3))); subplot(144) DATA = plt(FIELD(:,:,it4)); pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap gray xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it4))); % suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),... % ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']); save_figure end %% if 0 %% Show frame in kspace tf = 0; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600]) subplot(221); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$'); subplot(222); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) fftshift((abs(x)),2); FIELD = squeeze(Nipj(1,2,:,:,:)); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(FIELD(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{pj=01}|$'); subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$'); save_figure end %% if 0 %% Hermite energy spectra % tf = Ts2D(end-3); -fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1000, 300]); +fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]); plt = @(x) squeeze(x); for ij = 1:Nji subplotnum = 100+Nji*10+ij; subplot(subplotnum) for it5 = 1:2:Ns5D alpha = it5*1.0/Ns5D; loglog(Pi(1:2:end),plt(epsilon_i_pj(1:2:end,ij,it5)),... 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],... 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on; end grid on; xlabel('$p$'); TITLE = ['$\sum_{kr,kz} |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE); end save_figure end %% if 0 %% Laguerre energy spectra % tf = Ts2D(end-3); fig = figure; FIGNAME = ['laguerre_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]); plt = @(x) squeeze(x); for it5 = 1:2:Ns5D alpha = it5*1.0/Ns5D; loglog(Ji,plt(max(epsilon_i_pj(:,:,it5),[],1)),... 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],... 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on; end grid on; xlabel('$j$'); TITLE = ['$\max_p\sum_{kr,kz} |N_i^{pj}|^2$']; title(TITLE); save_figure end %% no_AA = (2:floor(2*Nkr/3)); -tKHI = 200; +tKHI = 100; [~,itKHI] = min(abs(Ts2D-tKHI)); after_KHI = (itKHI:Ns2D); if 0 %% Phi frequency space time diagram at kz=0 fig = figure; FIGNAME = ['phi_freq_diag_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]); [TY,TX] = meshgrid(Ts2D(after_KHI),kr(no_AA)); - pclr = pcolor(TX,TY,log10(squeeze(abs(PHI(no_AA,1,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar; + pclr = pcolor(TX,TY,(squeeze(abs(PHI(no_AA,1,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar; ylabel('$t c_s/R$'), xlabel('$0 0 % disp('Warning : dkperp < 0.0002'); % end %% %% We compute only on a kperp grid with dk space from 0 to kperpmax kperp = unique([0:dk:(sqrt(2)*kmax),sqrt(2)*kmax]); kperpmax = sqrt(2) * kmax; %% n_ = 1; for k_ = kperp + disp(['----------Computing matrix ',num2str(n_),'/',num2str(Nperp),'----------']) %% Script to run COSOlver in order to create needed collision matrices COSOlver_path = '../../Documents/MoliSolver/COSOlver/'; - COSOLVER.pmaxe = 10; - COSOLVER.jmaxe = 5; - COSOLVER.pmaxi = 10; - COSOLVER.jmaxi = 5; + COSOLVER.pmaxe = 20; + COSOLVER.jmaxe = 10; + COSOLVER.pmaxi = 20; + COSOLVER.jmaxi = 10; COSOLVER.kperp = k_; COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); % rule of thumb for sum truncation - COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2)); + COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); COSOLVER.idxT4max = 40; COSOLVER.neFLRs = 0; % ... only for GK abel COSOLVER.npeFLR = 0; % ... only for GK abel COSOLVER.niFLRs = 0; % ... only for GK abel COSOLVER.npiFLR = 0; % ... only for GK abel COSOLVER.gk = 1; COSOLVER.co = 3; % ! 0 = nothing % ! 1 = coulomb % ! 2 = pitch-angle (with constant coll.freq.) % ! 3 = sugama % ! 4 = pitch-angle with veloty dependent coll. freq. % ! 5 = improved sugama % ! 6 = Hirschman-Sigmar Clarke % ! 7 = Abel (only for like species) % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.) % ! 9 = GK coulomb polarization term write_fort90_COSOlver k_string = sprintf('%0.4f',k_); mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),... '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),... '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),... '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),... '_JE_12',... '_NFLR_',num2str(COSOLVER.neFLR),... '_kperp_',k_string,'.h5']; if (exist(mat_file_name,'file') > 0) disp(['Matrix available for kperp = ',k_string]); else cd ../../Documents/MoliSolver/COSOlver/ disp(['Matrix not found for kperp = ',k_string]); disp([num2str(n_),'/',Nperp]) disp('computing...'); CMD = 'mpirun -np 6 bin/CO 2 2 2 > out.txt'; disp(CMD); system(CMD); system(CMD); disp('..done'); cd ../../../HeLaZ/wk end + n_ = n_ + 1; end \ No newline at end of file diff --git a/wk/daint_run.m b/wk/daint_run.m index 772458f..e7f53ac 100644 --- a/wk/daint_run.m +++ b/wk/daint_run.m @@ -1,67 +1,89 @@ %clear all; addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss -CLUSTER.NODES = '1'; % number of nodes -CLUSTER.NTPN = '36'; % N tasks per node (mpi processes) -CLUSTER.NTPC = '1'; % N tasks per core (openmp threads) -CLUSTER.CPUPT = '1'; % CPU per task (number of CPU per mpi proc) +NP_P = 2; % MPI processes along p +NP_KR = 18; % MPI processes along kr CLUSTER.PART = 'normal'; % debug or normal if(strcmp(CLUSTER.PART,'debug')); CLUSTER.TIME = '00:30:00'; end; CLUSTER.MEM = '12GB'; % Memory CLUSTER.JNAME = 'gamma_inf';% Job name %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency -ETAB = 0.7; % Magnetic gradient -NU_HYP = 0.1; % Hyperdiffusivity coefficient +ETAB = 0.6; % Magnetic gradient +NU_HYP = 1.0; % Hyperdiffusivity coefficient %% GRID PARAMETERS N = 200; % Frequency gridpoints (Nkr = N/2) -L = 120; % Size of the squared frequency domain -P = 10; % Electron and Ion highest Hermite polynomial degree -J = 05; % Electron and Ion highest Laguerre polynomial degree +L = 120; % Size of the squared frequency domain +P = 12; % Electron and Ion highest Hermite polynomial degree +J = 06; % Electron and Ion highest Laguerre polynomial degree MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% TIME PARAMETERS -TMAX = 2500; % Maximal time unit +TMAX = 1000; % Maximal time unit DT = 1e-2; % Time step -SPS0D = 1; % Sampling per time unit for profiler -SPS2D = 1/2; % Sampling per time unit for 2D arrays -SPS5D = 1/10; % Sampling per time unit for 5D arrays -SPSCP = 0; % Sampling per time unit for checkpoints +SPS0D = 1; % Sampling per time unit for profiler +SPS2D = 1; % Sampling per time unit for 2D arrays +SPS5D = 1/40; % Sampling per time unit for 5D arrays +SPSCP = 0; % Sampling per time unit for checkpoints RESTART = 1; % To restart from last checkpoint -JOB2LOAD= 1; +JOB2LOAD= 2; %% OPTIONS -% SIMID = ['debug']; % Name of the simulation -SIMID = ['Daint_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation +SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation +% SIMID = 'test_marconi_sugama'; % Name of the simulation SIMID = sprintf(SIMID,NU); -CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty) +PREFIX =[]; +% PREFIX = sprintf('%d_%d_',NP_P, NP_KR); +% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) +CO = 1; CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P) +NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS) KERN = 0; % Kernel model (0 : GK) +INIT_PHI= 1; % Start simulation with a noisy phi and moments +%% OUTPUTS +W_DOUBLE = 1; +W_GAMMA = 1; +W_PHI = 1; +W_NA00 = 1; +W_NAPJ = 1; +W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% fixed parameters (for current study) KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; ETAT = 0.0; % Temperature gradient ETAN = 1.0; % Density gradient TAU = 1.0; % e/i temperature ratio +% Compute processes distribution +Ntot = NP_P * NP_KR; +Nnodes = ceil(Ntot/36); +Nppn = Ntot/Nnodes; +CLUSTER.NODES = num2str(Nnodes); % MPI process along p +CLUSTER.NTPN = num2str(Nppn); % MPI process along kr +CLUSTER.CPUPT = '1'; % CPU per task +CLUSTER.NTPC = '1'; % N tasks per core (openmp threads) %% Run file management scripts setup write_sbash_daint system('rm fort.90 setup_and_run.sh batch_script.sh'); disp('done'); +if(mod(NP_P*NP_KR,36)~= 0) + disp('WARNING : unused cores (ntot cores must be a 36 multiple)'); +end \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index a529cd2..6812cba 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,62 +1,64 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 1.0; % Collision frequency -ETAB = 0.5; % Magnetic gradient +ETAB = 0.6; % Magnetic gradient ETAN = 1.0; % Density gradient NU_HYP = 1.0; %% GRID PARAMETERS N = 200; % Frequency gridpoints (Nkr = N/2) -L = 120; % Size of the squared frequency domain -PMAXE = 10; % Highest electron Hermite polynomial degree -JMAXE = 05; % Highest '' Laguerre '' -PMAXI = 10; % Highest ion Hermite polynomial degree -JMAXI = 05; % Highest '' Laguerre '' +L = 80; % Size of the squared frequency domain +PMAXE = 02; % Highest electron Hermite polynomial degree +JMAXE = 01; % Highest '' Laguerre '' +PMAXI = 02; % Highest ion Hermite polynomial degree +JMAXI = 01; % Highest '' Laguerre '' +MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f +MU_J = 0.0/JMAXI^2; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% TIME PARAMETERS -TMAX = 50; % Maximal time unit +TMAX = 2500; % Maximal time unit DT = 2e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1/2; % Sampling per time unit for 2D arrays SPS5D = 1/4; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints/10 -RESTART = 0; % To restart from last checkpoint -JOB2LOAD= 0; +RESTART = 1; % To restart from last checkpoint +JOB2LOAD= 3; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) -CO = 2 ; +CO = 1; CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P) NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) -SIMID = 'test_SGGK_full_size'; % Name of the simulation -NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) +SIMID = 'test_stability'; % Name of the simulation +% SIMID = ['local_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation +% SIMID = sprintf(SIMID,NU); +NON_LIN = 1; % activate non-linearity (is cancelled if KREQ0 = 1) %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused KERN = 0; % Kernel model (0 : GK) KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; -kmax = N*pi/L;% Highest fourier mode +kmax = 2/3*N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio % kmaxcut = 2.5; MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient -MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f -MU_J = 0.0/JMAXI^3; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f INIT_PHI= 1; % Start simulation with a noisy phi and moments %% Setup and file management setup system('rm fort.90'); \ No newline at end of file diff --git a/wk/marconi_run.m b/wk/marconi_run.m index 450f2d1..346e96e 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,86 +1,87 @@ %clear all; addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS -CLUSTER.TIME = '00:10:00'; % allocation time hh:mm:ss -CLUSTER.PART = 'dbg'; % dbg or prod -CLUSTER.MEM = '16GB'; % Memory +CLUSTER.PART = 'prod'; % dbg or prod +CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss +if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end; +CLUSTER.MEM = '64GB'; % Memory CLUSTER.JNAME = 'gamma_inf';% Job name -NP_P = 1; % MPI processes along p -NP_KR = 1; % MPI processes along kr +NP_P = 2; % MPI processes along p +NP_KR = 24; % MPI processes along kr %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency ETAB = 0.6; % Magnetic gradient NU_HYP = 1.0; % Hyperdiffusivity coefficient %% GRID PARAMETERS N = 200; % Frequency gridpoints (Nkr = N/2) -L = 120; % Size of the squared frequency domain -P = 04; % Electron and Ion highest Hermite polynomial degree -J = 04; % Electron and Ion highest Laguerre polynomial degree +L = 80; % Size of the squared frequency domain +P = 12; % Electron and Ion highest Hermite polynomial degree +J = 06; % Electron and Ion highest Laguerre polynomial degree MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% TIME PARAMETERS -TMAX = 120; % Maximal time unit -DT = 2e-2; % Time step +TMAX = 4000; % Maximal time unit +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler -SPS2D = 1; % Sampling per time unit for 2D arrays -SPS5D = 1/40; % Sampling per time unit for 5D arrays +SPS2D = 1/2; % Sampling per time unit for 2D arrays +SPS5D = 1/100; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints -RESTART = 0; % To restart from last checkpoint -JOB2LOAD= 0; +RESTART = 1; % To restart from last checkpoint +JOB2LOAD= 4; %% OPTIONS -% SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation -SIMID = 'test_marconi_sugama'; % Name of the simulation +SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation SIMID = sprintf(SIMID,NU); +% SIMID = 'test_marconi_sugama'; % Name of the simulation PREFIX =[]; % PREFIX = sprintf('%d_%d_',NP_P, NP_KR); % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) -CO = 2; +CO = 1; CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P) NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_PHI= 1; % Start simulation with a noisy phi and moments %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% fixed parameters (for current study) KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; ETAT = 0.0; % Temperature gradient ETAN = 1.0; % Density gradient TAU = 1.0; % e/i temperature ratio % Compute processes distribution Ntot = NP_P * NP_KR; Nnodes = ceil(Ntot/48); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kr CLUSTER.CPUPT = '1'; % CPU per task %% Run file management scripts setup write_sbash_marconi system('rm fort.90 setup_and_run.sh batch_script.sh'); disp('done'); if(mod(NP_P*NP_KR,48)~= 0) disp('WARNING : unused cores (ntot cores must be a 48 multiple)'); end \ No newline at end of file