diff --git a/matlab/create_gif_1D.m b/matlab/create_gif_1D.m index 5f7369f..296117f 100644 --- a/matlab/create_gif_1D.m +++ b/matlab/create_gif_1D.m @@ -1,50 +1,52 @@ 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))); flag = 0; if hmax == hmin disp('Warning : h = hmin = hmax = const') else % Setup figure frame fig = figure; 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)); + scale = max(FIELD(:,n))*SCALING + (1-SCALING); 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; end disp(' ') disp(['Gif saved @ : ',GIFNAME]) end diff --git a/matlab/setup.m b/matlab/setup.m index d9beeda..e5ac181 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.Nr = N; % r grid resolution GRID.Lr = L; % r length GRID.Nz = N * (1-KREQ0) + KREQ0; % z '' GRID.Lz = L * (1-KREQ0); % z '' GRID.kpar = KPAR; % Model parameters MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end; MODEL.mu = MU; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % mass ratio sqrt(m_a/m_i) MODEL.sigma_e = 0.0233380; MODEL.sigma_i = 1.0; % charge q_a/e MODEL.q_e =-1.0; MODEL.q_i = 1.0; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.eta_n = ETAN; % source term kappa for HW MODEL.eta_T = ETAT; % Temperature MODEL.eta_B = ETAB; % Magnetic MODEL.lambdaD = LAMBDAD; % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; if INIT_PHI; INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end; INITIAL.initback_moments = 0.0e-5; INITIAL.initnoise_moments = NOISE0; INITIAL.iseed = 42; INITIAL.selfmat_file = '''null'''; INITIAL.eimat_file = '''null'''; INITIAL.iemat_file = '''null'''; if (CO == -3) % Write matrice filename for Full Coulomb DK cmat_pmaxe = 25; cmat_jmaxe = 18; cmat_pmaxi = 25; cmat_jmaxi = 18; INITIAL.selfmat_file = ... ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_pamaxx_10.h5''']; INITIAL.eimat_file = ... ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5''']; INITIAL.iemat_file = ... ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5''']; elseif (CO == -2) % Write matrice filename for DK Sugama cmat_pmaxe = 10; cmat_jmaxe = 5; cmat_pmaxi = 10; cmat_jmaxi = 5; INITIAL.selfmat_file = ... ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_JE_12.h5''']; INITIAL.eimat_file = ... ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5''']; INITIAL.iemat_file = ... ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5''']; elseif (CO == 2) % Write matrice filename for Sugama GK cmat_pmaxe = 10; cmat_jmaxe = 5; cmat_pmaxi = 10; cmat_jmaxi = 5; INITIAL.selfmat_file = ... ['''../../../iCa/self_Coll_GKE_1_GKI_1_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... - num2str(cmat_jmaxi),'_JE_12_NFLR_5_''']; + num2str(cmat_jmaxi),'_JE_12_''']; INITIAL.eimat_file = ... ['''../../../iCa/ei_Coll_GKE_1_GKI_1_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... - num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_NFLRe_5_NFLRi_5_''']; + num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_''']; INITIAL.iemat_file = ... ['''../../../iCa/ie_Coll_GKE_1_GKI_1_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),... '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',... - num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_NFLRe_5_NFLRi_5_''']; + num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_''']; elseif (CO == 3) % Full Coulomb GK disp('Warning, FCGK not implemented yet') elseif (CO == -1) % DGDK disp('Warning, DGDK not debugged') end % Naming and creating input file if (CO == -3); CONAME = 'FCDK'; elseif(CO == -2); CONAME = 'SGDK'; elseif(CO == -1); CONAME = 'DGDK'; elseif(CO == 0); CONAME = 'LB'; elseif(CO == 1); CONAME = 'DGGK'; elseif(CO == 2); CONAME = 'SGGK'; elseif(CO == 3); CONAME = 'FCGK'; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end if (PMAXE == PMAXI) && (JMAXE == JMAXI) degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; else degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; end degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',... CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_']; gridname = ['L_',num2str(L),'_']; if (exist('PREFIX','var') == 0); PREFIX = []; end; if (exist('SUFFIX','var') == 0); SUFFIX = []; end; PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end; OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT); if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; OUTPUTS.resfile0 = '''outputs'''; OUTPUTS.rstfile0 = '''checkpoint'''; OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) end if ~exist(BASIC.RESDIR, 'dir') mkdir(BASIC.RESDIR) end %% Compile and WRITE input file INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; system(MAKE); %% disp(['Set up ',SIMID]); disp([resolution,gridname,degngrad]); if RESTART disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/matlab/write_fort90_COSOlver.m b/matlab/write_fort90_COSOlver.m index 46da401..066f825 100644 --- a/matlab/write_fort90_COSOlver.m +++ b/matlab/write_fort90_COSOlver.m @@ -1,84 +1,84 @@ INPUT = [COSOlver_path,'fort.90']; fid = fopen(INPUT,'wt'); fprintf(fid,'&BASIC\n'); fprintf(fid,['pmaxe = ', num2str(COSOLVER.pmaxe),'\n']); fprintf(fid,['jmaxe = ', num2str(COSOLVER.jmaxe),'\n']); fprintf(fid,['pmaxi = ', num2str(COSOLVER.pmaxi),'\n']); fprintf(fid,['jmaxi = ', num2str(COSOLVER.jmaxi),'\n']); fprintf(fid,'JEmaxx=12\n'); fprintf(fid,'PMmaxx=16\n'); fprintf(fid,'\n'); fprintf(fid,['neFLR=',num2str(COSOLVER.neFLR),'\n']); fprintf(fid,['neFLRs =',num2str(COSOLVER.neFLRs),'\n']); fprintf(fid,['npeFLR=',num2str(COSOLVER.npeFLR),'\n']); fprintf(fid,['niFLR=',num2str(COSOLVER.niFLR),'\n']); fprintf(fid,['niFLRs=',num2str(COSOLVER.niFLRs),'\n']); fprintf(fid,['npiFLR=',num2str(COSOLVER.npiFLR),'\n']); fprintf(fid,'\n'); fprintf(fid,['eecolls=.true.','\n']); fprintf(fid,['iicolls=.true.','\n']); fprintf(fid,['eicolls=.true.','\n']); fprintf(fid,['iecolls=.true.','\n']); fprintf(fid,'/\n'); fprintf(fid,'\n'); fprintf(fid,'&BASIS_TRANSFORMATION_PAR\n'); fprintf(fid,['T5dir = ','''/misc/coeffs_backup/T5src/''','\n']); fprintf(fid,['T4dir = ','''/misc/T4/NNT4_L000x200_K000x200_P000x200_J000x200/''','\n']); -fprintf(fid,'idxT4max = 30\n'); +fprintf(fid,['idxT4max = ',num2str(COSOLVER.idxT4max),'\n']); fprintf(fid,'idxT5max = 0\n'); fprintf(fid,'IFT4 = .true.\n'); fprintf(fid,'IFT5 = .false.\n'); fprintf(fid,'/\n'); fprintf(fid,'\n'); fprintf(fid,'&MODEL_PAR \n'); fprintf(fid,'nu=1\n'); fprintf(fid,'mu=0.023338\n'); fprintf(fid,'tau=1\n'); fprintf(fid,['kperp=',num2str(COSOLVER.kperp,16),'\n']); fprintf(fid,'/\n'); fprintf(fid,'\n'); fprintf(fid,'&OPERATOR_MODEL \n'); fprintf(fid,['ETEST=', num2str(COSOLVER.co),'\n']); fprintf(fid,['EBACK=', num2str(COSOLVER.co),'\n']); fprintf(fid,['ITEST=', num2str(COSOLVER.co),'\n']); fprintf(fid,['IBACK=', num2str(COSOLVER.co),'\n']); fprintf(fid,'\n'); fprintf(fid,['ESELF=', num2str(COSOLVER.co),'\n']); fprintf(fid,['ISELF=', num2str(COSOLVER.co),'\n']); fprintf(fid,'\n'); fprintf(fid,['GKE = ', num2str(COSOLVER.gk),'\n']); fprintf(fid,['GKI = ', num2str(COSOLVER.gk),'\n']); fprintf(fid,'DKTEST = .F.\n'); fprintf(fid,'DKBACK = .F.\n'); fprintf(fid,'ADDTEST = .T.\n'); fprintf(fid,'ADDBACK = .T.\n'); fprintf(fid,'only_gk_part = .false.\n'); fprintf(fid,'only_symmetric = .false.\n'); fprintf(fid,'/\n'); fprintf(fid,'\n'); fprintf(fid,'&MPI\n'); fprintf(fid,'Nprocs_j_ii = 1\n'); fprintf(fid,'MPI_balanced = .false.\n'); fprintf(fid,'/\n'); fprintf(fid,'\n'); fprintf(fid,'&OUTPUT_PAR\n'); fprintf(fid,'OVERWRITE=.true.\n'); fprintf(fid,'nsave_ei=0\n'); fprintf(fid,'nsave_ie=0\n'); fprintf(fid,'nsave_ee=0\n'); fprintf(fid,'nsave_ii=2\n'); fprintf(fid,'ifrestart=.false.\n'); fprintf(fid,['suffix_resfile=','''''','\n']); fprintf(fid,['outdir = ','''../../../HeLaZ/iCa''','\n']); fprintf(fid,'/\n'); fprintf(fid,'\n'); fclose(fid); \ No newline at end of file diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index 6d43bb3..50769e0 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,415 +1,454 @@ %% Load results outfile =''; -if 0 +if 1 %% Load from Marconi outfile =''; -outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_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.4_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.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_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'; BASIC.RESDIR = load_marconi(outfile); end if 0 %% Load from Daint outfile =''; 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 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)))); 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); 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(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}$') subplot(223) plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on; grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show'); 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)]); +suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]); 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); -time_array = [1, 100, 400, 1000]; fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1000, 300]); 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 |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE); + 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; +[~,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; + ylabel('$t c_s/R$'), xlabel('$0 0 - disp('Warning : dkperp < 0.0002'); -end +% %% Check if the differences btw kperp is larger than naming precision +% dkperp = diff(kperp); +% warning = sum(dkperp<0.0002); +% if warning > 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 %% 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.kperp = k_; - COSOLVER.neFLR = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation - COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2)); + 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.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([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 end \ No newline at end of file diff --git a/wk/linear_study.m b/wk/linear_study.m index 9e3b3b6..9734c39 100644 --- a/wk/linear_study.m +++ b/wk/linear_study.m @@ -1,209 +1,216 @@ +for NU = [1.0 0.1 0.01] +for ETAB = [0.5 0.6 0.7 0.8] +for CO = [-3 -2 -1 0 1 2] %clear all; addpath(genpath('../matlab')) % ... add default_plots_options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.1; % Collision frequency +% NU = 1.0; % Collision frequency TAU = 1.0; % e/i temperature ratio -ETAB = 0.6; +% ETAB = 0.5; ETAN = 1.0; % Density gradient ETAT = 0.0; % Temperature gradient NU_HYP = 0.0; % Hyperdiffusivity coefficient LAMBDAD = 0.0; NOISE0 = 1.0e-5; %% GRID PARAMETERS -N = 10; % Frequency gridpoints (Nkr = N/2) +N = 100; % Frequency gridpoints (Nkr = N/2) L = 120; % Size of the squared frequency domain KREQ0 = 1; % put kr = 0 MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% TIME PARMETERS TMAX = 200; % Maximal time unit -DT = 1e-2; % Time step +DT = 2e-2; % Time step SPS0D = 0.5; % Sampling per time unit for 2D arrays SPS2D = 1; % Sampling per time unit for 2D arrays SPS5D = 1/2; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints RESTART = 0; % To restart from last checkpoint JOB2LOAD= 00; %% OPTIONS -SIMID = 'linear_study_SugamaGK'; % Name of the simulation +SIMID = 'linear_study'; % Name of the simulation NON_LIN = 0 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) -CO = 2; +% CO = 2; CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P) NL_CLOS = 0; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_PHI= 0; % Start simulation with a noisy phi and moments %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused % DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1) -JOBNUM = 00; +JOBNUM = 00; KPAR = 0.0; % Parellel wave vector component HD_CO = 0.5; % Hyper diffusivity cutoff ratio kmax = N*pi/L;% Highest fourier mode MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient %% PARAMETER SCANS if 1 %% Parameter scan over PJ PA = [2, 3, 4, 6, 8, 10]; JA = [1, 2, 2, 3, 4, 5]; DTA= DT./sqrt(JA); mup_ = MU_P; muj_ = MU_J; % PA = [4]; % JA = [1]; Nparam = numel(PA); param_name = 'PJ'; gamma_Ni00 = zeros(Nparam,floor(N/2)+1); -gamma_Ni21 = zeros(Nparam,floor(N/2)+1); +gamma_Nipj = zeros(Nparam,floor(N/2)+1); Bohm_transport = zeros(Nparam,1); Ni00_ST = zeros(Nparam,floor(N/2)+1,SPS2D*TMAX); for i = 1:Nparam % Change scan parameter PMAXE = PA(i); PMAXI = PA(i); JMAXE = JA(i); JMAXI = JA(i); DT = DTA(i); MU_P = mup_/PMAXE^2; MU_J = muj_/JMAXE^3; setup % Run linear simulation % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz 1 1; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 2 3; cd ../../../wk']) % Load and process results load_results tend = Ts2D(end); tstart = 0.4*tend; for ikr = 1:N/2+1 gamma_Ni00(i,ikr) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,1,:))),tstart,tend); Ni00_ST(i,ikr,1:numel(Ts2D)) = squeeze((Ni00(ikr,1,:))); end tend = Ts5D(end); tstart = 0.4*tend; for ikr = 1:N/2+1 - gamma_Ni21(i,ikr) = LinearFit_s(Ts5D,squeeze(abs(Nipj(3,2,ikr,1,:))),tstart,tend); + gamma_Nipj(i,ikr) = LinearFit_s(Ts5D,squeeze(max(max(abs(Nipj(:,:,ikr,1,:)),[],1),[],2)),tstart,tend); end gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0)); - gamma_Ni21(i,:) = real(gamma_Ni21(i,:) .* (gamma_Ni21(i,:)>=0.0)); + gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0)); % kzmax = abs(kr(ikzmax)); % Bohm_transport(i) = ETAB/ETAN*gmax/kzmax^2; % Clean output system(['rm -r ',BASIC.RESDIR]) end if 1 %% Plot SCALE = 1;%sqrt(2); fig = figure; FIGNAME = 'linear_study'; plt = @(x) x; subplot(211) for i = 1:Nparam clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); plot(plt(SCALE*kr),plt(gamma_Ni00(i,:)),... 'Color',clr,... 'LineStyle',linestyle{1},... 'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]); hold on; end grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})\rho_s/c_s$'); xlim([0.0,max(kr)]); title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME]) legend('show') subplot(212) for i = 1:Nparam clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); - plot(plt(SCALE*kr),plt(gamma_Ni21(i,:)),... + plot(plt(SCALE*kr),plt(gamma_Nipj(i,:)),... 'Color',clr,... 'LineStyle',linestyle{1},... 'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]); hold on; end - grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{21})\rho_s/c_s$'); xlim([0.0,max(kr)]); + grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(\max_{pj}N_i^{pj})\rho_s/c_s$'); xlim([0.0,max(kr)]); title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME]) legend('show') saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']); saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']); end end if 0 %% Parameter scan over eta_B PMAXE = 6; PMAXI = 6; JMAXE = 3; JMAXI = 3; TMAX = 400; DT = 0.001; eta_B = [0.45, 0.5, 0.55, 0.6, 0.65]; Nparam = numel(eta_B); param_name = 'etaB'; CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) NU = 1e-2; % Collision frequency Bohm_transport = zeros(Nparam,1); gamma_Ni = zeros(Nparam,N); for i = 1:Nparam % Change scan parameter ETAB = eta_B(i); setup % Run linear simulation system(... ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ./../../../bin/helaz; cd ../../../wk']... ) % Load and process results load_results tend = Ts2D(end); tstart = 0.4*tend; for ikr = 1:N/2+1 gamma_Ni(i,ikr) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,1,:))),tstart,tend); Ni00_ST(i,ikr,1:numel(Ts2D)) = squeeze((Ni00(ikr,1,:))); end gamma_Ni(i,:) = real(gamma_Ni(i,:) .* (gamma_Ni(i,:)>=0.0)); [gmax,ikzmax] = max(gamma_Ni(i,:)); kzmax = abs(kr(ikzmax)); Bohm_transport(i) = ETAB/ETAN*gmax/kzmax^2; % Clean output system(['rm -r ',BASIC.RESDIR]) end if 0 %% Plot fig = figure; FIGNAME = 'linear_study'; plt = @(x) circshift(x,N/2-1); for i = 1:Nparam clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); plot(plt(kz),plt(gamma_Ni(i,:)),... 'Color',clr,... 'LineStyle',linestyle{1},... 'DisplayName',['$\eta_B=$',num2str(eta_B(i))]); hold on; end grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma(N_i^{00})\rho_2/c_s$'); xlim([0.0,max(kz)]); title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) legend('show') saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']); end if 1 %% Plot fig = figure; FIGNAME = 'mixing_length'; plot(eta_B, Bohm_transport) grid on; xlabel('$L_n/L_B$'); ylabel('$\eta\gamma_{max}/k_{max}^2$'); title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']); end +%% end +end +end +end \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index aa82246..a529cd2 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,65 +1,62 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.01; % Collision frequency -ETAB = 0.8; % Magnetic gradient +NU = 1.0; % Collision frequency +ETAB = 0.5; % Magnetic gradient ETAN = 1.0; % Density gradient NU_HYP = 1.0; %% GRID PARAMETERS -N = 50; % Frequency gridpoints (Nkr = N/2) -L = 100; % Size of the squared frequency domain -PMAXE = 4; % Highest electron Hermite polynomial degree -JMAXE = 4; % Highest '' Laguerre '' -PMAXI = 4; % Highest ion Hermite polynomial degree -JMAXI = 4; % Highest '' Laguerre '' +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 '' %% TIME PARAMETERS -TMAX = 2000; % Maximal time unit +TMAX = 50; % 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; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) -CO = 1; +CO = 2 ; 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) -KERN = 0; % Kernel model (0 : GK) -INIT_PHI= 1; % Start simulation with a noisy phi and moments -% SIMID = ['local_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation -% SIMID = sprintf(SIMID,NU); -% SIMID = 'test_init_phi'; % Name of the simulation -SIMID = ['test_nonlin_NL_CLOS_',num2str(NL_CLOS)]; % Name of the simulation +SIMID = 'test_SGGK_full_size'; % Name of the simulation +NON_LIN = 1 *(1-KREQ0); % 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; -NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) kmax = 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 0d5ab39..450f2d1 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,86 +1,86 @@ %clear all; addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS -CLUSTER.TIME = '12:00:00'; % allocation time hh:mm:ss -CLUSTER.PART = 'prod'; % dbg or prod +CLUSTER.TIME = '00:10:00'; % allocation time hh:mm:ss +CLUSTER.PART = 'dbg'; % dbg or prod CLUSTER.MEM = '16GB'; % Memory CLUSTER.JNAME = 'gamma_inf';% Job name -NP_P = 2; % MPI processes along p -NP_KR = 24; % MPI processes along kr +NP_P = 1; % MPI processes along p +NP_KR = 1; % 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 = 10; % Electron and Ion highest Hermite polynomial degree -J = 05; % Electron and Ion highest Laguerre polynomial degree +P = 04; % Electron and Ion highest Hermite polynomial degree +J = 04; % 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 = 250; % Maximal time unit -DT = 5e-4; % Time step -SPS0D = 1; % Sampling per time unit for profiler -SPS2D = 1; % Sampling per time unit for 2D arrays -SPS5D = 1/50; % Sampling per time unit for 5D arrays -SPSCP = 0; % Sampling per time unit for checkpoints -RESTART = 1; % To restart from last checkpoint -JOB2LOAD= 1; +TMAX = 120; % Maximal time unit +DT = 2e-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 +SPSCP = 0; % Sampling per time unit for checkpoints +RESTART = 0; % To restart from last checkpoint +JOB2LOAD= 0; %% OPTIONS -SIMID = ['HeLaZ_v2.4_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation -% SIMID = 'Marconi_parallel_scaling_2D'; % 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); PREFIX =[]; % PREFIX = sprintf('%d_%d_',NP_P, NP_KR); % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK) -CO = 1; +CO = 2; 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 24 multiple)'); + disp('WARNING : unused cores (ntot cores must be a 48 multiple)'); end \ No newline at end of file