diff --git a/matlab/compile_results.m b/matlab/compile_results.m index b4c3bae..dc515ac 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,159 +1,212 @@ +function [DATA] = compile_results(DIRECTORY,JOBNUMMIN,JOBNUMMAX) +DATA = {}; CONTINUE = 1; JOBNUM = JOBNUMMIN; JOBFOUND = 0; -TJOB_SE = []; % Start and end times of jobs -NU_EVOL = []; % evolution of parameter nu between jobs -CO_EVOL = []; % evolution of CO -MU_EVOL = []; % evolution of parameter mu between jobs -K_N_EVOL = []; % -L_EVOL = []; % -DT_EVOL = []; % +DATA.TJOB_SE = []; % Start and end times of jobs +DATA.NU_EVOL = []; % evolution of parameter nu between jobs +DATA.CO_EVOL = []; % evolution of CO +DATA.MU_EVOL = []; % evolution of parameter mu between jobs +DATA.K_N_EVOL = []; % +DATA.L_EVOL = []; % +DATA.DT_EVOL = []; % % FIELDS Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; HFLUX_ = []; GGAMMA_ = []; PGAMMA_ = []; PHI_ = []; DENS_E_ = []; DENS_I_ = []; TEMP_E_ = []; TEMP_I_ = []; Ts0D_ = []; Ts3D_ = []; Ts5D_ = []; Sipj_ = []; Sepj_ = []; Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old; Pi_max=0; Pe_max=0; Ji_max=0; Je_max=0; while(CONTINUE) - filename = sprintf([BASIC.MISCDIR,'outputs_%.2d.h5'],JOBNUM); + filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM); if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX) % Load results of simulation #JOBNUM - load_results +% load_results + %% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + disp(['Loading ',filename]) + % Loading from output file + CPUTIME = h5readatt(filename,'/data/input','cpu_time'); + DT_SIM = h5readatt(filename,'/data/input','dt'); + [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename); + W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y'); + W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y'); + W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y'); + W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y'); + W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y'); + W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y'); + W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y'); + W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y'); +% KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y'); + KIN_E = 1; + % Check polynomials degrees Pe_new= numel(Pe); Je_new= numel(Je); Pi_new= numel(Pi); Ji_new= numel(Ji); if(Pe_max < Pe_new); Pe_max = Pe_new; end; if(Je_max < Je_new); Je_max = Je_new; end; if(Pi_max < Pi_new); Pi_max = Pi_new; end; if(Ji_max < Ji_new); Ji_max = Ji_new; end; % If a degree is larger than previous job, put them in a larger array if (sum([Pe_new, Je_new, Pi_new, Ji_new]>[Pe_old, Je_old, Pi_old, Ji_old]) >= 1) if W_NAPJ tmp = Nipj_; sz = size(tmp); Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; tmp = Nepj_; sz = size(tmp); Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; end if W_SAPJ tmp = Sipj_; sz = size(tmp); Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; tmp = Sepj_; sz = size(tmp); Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; end % If a degree is smaller than previous job, put zero to add. deg. elseif (sum([Pe_new, Je_new, Pi_new, Ji_new]<[Pe_old, Je_old, Pi_old, Ji_old]) >= 1 && Pe_old ~= 1e9) if W_NAPJ tmp = Nipj; sz = size(tmp); Nipj = zeros(cat(1,[Pi_max,Ji_max]',sz(3:end)')'); Nipj(1:Pi_new,1:Ji_new,:,:,:) = tmp; tmp = Nepj; sz = size(tmp); Nepj = zeros(cat(1,[Pe_max,Je_max]',sz(3:end)')'); Nepj(1:Pe_new,1:Je_new,:,:,:) = tmp; end if W_SAPJ tmp = Sipj; sz = size(tmp); Sipj = zeros(cat(1,[Pi_max,Ji_max]',sz(3:end)')'); Sipj(1:Pi_new,1:Ji_new,:,:,:) = tmp; tmp = Sepj; sz = size(tmp); Sepj = zeros(cat(1,[Pe_max,Je_max]',sz(3:end)')'); Sepj(1:Pe_new,1:Je_new,:,:,:) = tmp; end end - if W_GAMMA || W_HF - Ts0D_ = cat(1,Ts0D_,Ts0D); - end if W_GAMMA - GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); - PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); + [ GGAMMA_RI, Ts0D, ~] = load_0D_data(filename, 'gflux_ri'); + PGAMMA_RI = load_0D_data(filename, 'pflux_ri'); + GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI); clear GGAMMA_RI + PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI); clear PGAMMA_RI end if W_HF - HFLUX_ = cat(1,HFLUX_,HFLUX_X); + [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x'); + HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X end - if W_PHI || W_NA00 - Ts3D_ = cat(1,Ts3D_,Ts3D); - end if W_PHI - PHI_ = cat(4,PHI_,PHI); + [ PHI, Ts3D, ~] = load_3D_data(filename, 'phi'); + PHI_ = cat(4,PHI_,PHI); clear PHI end if W_NA00 if KIN_E - Ne00_ = cat(4,Ne00_,Ne00); + Ne00 = load_3D_data(filename, 'Ne00'); + Ne00_ = cat(4,Ne00_,Ne00); clear Ne00 end - Ni00_ = cat(4,Ni00_,Ni00); + [Ni00, Ts3D, ~] = load_3D_data(filename, 'Ni00'); + Ni00_ = cat(4,Ni00_,Ni00); clear Ni00 end if W_DENS if KIN_E - DENS_E_ = cat(4,DENS_E_,DENS_E); + [DENS_E, Ts3D, ~] = load_3D_data(filename, 'dens_e'); + DENS_E_ = cat(4,DENS_E_,DENS_E); clear DENS_E end - DENS_I_ = cat(4,DENS_I_,DENS_I); + [DENS_I, Ts3D, ~] = load_3D_data(filename, 'dens_i'); + DENS_I_ = cat(4,DENS_I_,DENS_I); clear DENS_I end if W_TEMP if KIN_E - TEMP_E_ = cat(4,TEMP_E_,TEMP_E); + [TEMP_E, Ts3D, ~] = load_3D_data(filename, 'temp_e'); + TEMP_E_ = cat(4,TEMP_E_,TEMP_E); clear TEMP_E end - TEMP_I_ = cat(4,TEMP_I_,TEMP_I); - end - if W_NAPJ || W_SAPJ - Ts5D_ = cat(1,Ts5D_,Ts5D); + [TEMP_I, Ts3D, ~] = load_3D_data(filename, 'temp_i'); + TEMP_I_ = cat(4,TEMP_I_,TEMP_I); clear TEMP_I end + if W_NAPJ - Nipj_ = cat(6,Nipj_,Nipj); + [Nipj, Ts5D, ~] = load_5D_data(filename, 'moments_i'); + Nipj_ = cat(6,Nipj_,Nipj); clear Nipj if KIN_E - Nepj_ = cat(6,Nepj_,Nepj); + Nepj = load_5D_data(filename, 'moments_e'); + Nepj_ = cat(6,Nepj_,Nepj); clear Nepj end end if W_SAPJ Sipj_ = cat(6,Sipj_,Sipj); if KIN_E Sepj_ = cat(6,Sepj_,Sepj); end end + Ts0D_ = cat(1,Ts0D_,Ts0D); + Ts3D_ = cat(1,Ts3D_,Ts3D); + Ts5D_ = cat(1,Ts5D_,Ts5D); % Evolution of simulation parameters load_params - TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)]; - NU_EVOL = [NU_EVOL NU NU]; - CO_EVOL = [CO_EVOL CO CO]; - MU_EVOL = [MU_EVOL MU MU]; - K_N_EVOL = [K_N_EVOL K_N K_N]; - L_EVOL = [L_EVOL L L]; - DT_EVOL = [DT_EVOL DT_SIM DT_SIM]; + DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)]; + DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU]; + DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO]; + DATA.MU_EVOL = [DATA.MU_EVOL DATA.MU DATA.MU]; + DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N]; + DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L]; + DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM]; JOBFOUND = JOBFOUND + 1; LASTJOB = JOBNUM; Pe_old = Pe_new; Je_old = Je_new; Pi_old = Pi_new; Ji_old = Ji_new; elseif (JOBNUM > JOBNUMMAX) CONTINUE = 0; disp(['found ',num2str(JOBFOUND),' results']); end JOBNUM = JOBNUM + 1; end -GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; HFLUX_X = HFLUX_; Ts0D = Ts0D_; -Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_; -Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts3D = Ts3D_; -DENS_E = DENS_E_; DENS_I = DENS_I_; TEMP_E = TEMP_E_; TEMP_I = TEMP_I_; -clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_ HFLUX_ -Sipj = Sipj_; Sepj = Sepj_; -clear Sipj_ Sepj_ +%% Build grids +Nkx = numel(kx); Nky = numel(ky); +[KY,KX] = meshgrid(ky,kx); +Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky); +dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1); +KPERP2 = KY.^2+KX.^2; +[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky)); + +Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z); +Lx = 2*pi/dkx; Ly = 2*pi/dky; +dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz; +x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x); +y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y); +%% Add everything in output structure +% Fields +DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_; +DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_; +if(KIN_E) +DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_; +end +DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_; +DATA.PHI = PHI_; +% grids +DATA.Pe = Pe; DATA.Pi = Pi; +DATA.kx = kx; DATA.ky = ky; DATA.z = z; +DATA.x = x; DATA.y = y; +DATA.ikx0 = ikx0; DATA.iky0 = iky0; +DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky; +DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji); +DATA.dir = DIRECTORY; +DATA.localdir = ['..',DIRECTORY(20:end)]; + JOBNUM = LASTJOB; -filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); + +filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM); +end \ No newline at end of file diff --git a/matlab/create_gif.m b/matlab/create_film.m similarity index 52% rename from matlab/create_gif.m rename to matlab/create_film.m index 1d0d571..b48cf20 100644 --- a/matlab/create_gif.m +++ b/matlab/create_film.m @@ -1,76 +1,103 @@ -title1 = GIFNAME; -FIGDIR = BASIC.RESDIR; -GIFNAME = [FIGDIR, GIFNAME,'.gif']; - -XNAME = latexize(XNAME); -YNAME = latexize(YNAME); -FIELDNAME = latexize(FIELDNAME); +function create_film(DATA,OPTIONS,format) +%% Plot options +FPS = 30; DELAY = 1/FPS; +INTERP = OPTIONS.INTERP; BWR = 1; NORMALIZED = 1; +T = DATA.Ts3D; +%% Processing +toplot = process_field(DATA,OPTIONS); +%% +FILENAME = [toplot.FILENAME,format]; +XNAME = ['$',toplot.XNAME,'$']; +YNAME = ['$',toplot.YNAME,'$']; +FIELDNAME = ['$',toplot.FIELDNAME,'$']; +FIELD = toplot.FIELD; X = toplot.X; Y = toplot.Y; +FRAMES = toplot.FRAMES; +switch format + case '.avi' + vidfile = VideoWriter(FILENAME,'Uncompressed AVI'); + vidfile.FrameRate = FPS; + open(vidfile); +end % Set colormap boundaries hmax = max(max(max(FIELD(:,:,FRAMES)))); hmin = min(min(min(FIELD(:,:,FRAMES)))); scale = -1; flag = 0; if hmax == hmin disp('Warning : h = hmin = hmax = const') else % Setup figure frame -fig = figure('Color','white','Position', [100, 100, 400, 400]); +fig = figure('Color','white','Position', toplot.DIMENSIONS); pcolor(X,Y,FIELD(:,:,1)); % to set up if BWR colormap(bluewhitered) end axis tight manual % this ensures that getframe() returns a consistent size if INTERP shading interp; end in = 1; nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); for n = FRAMES % loop over selected frames scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize if NORMALIZED pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\ caxis([-1,1]); else pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\ if CONST_CMAP caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map else caxis([-1,1]*scale); % adaptive color map end title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... ,', scale = ',sprintf('%.1e',scale)]); end title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... ,', scaling = ',sprintf('%.1e',scale)]); if INTERP shading interp; end - set(pclr, 'edgecolor','none'); axis square; + set(pclr, 'edgecolor','none'); pbaspect(toplot.ASPECT); if BWR colormap(bluewhitered) end xlabel(XNAME); ylabel(YNAME); %colorbar; drawnow % Capture the plot as an image frame = getframe(fig); - im = frame2im(frame); - [imind,cm] = rgb2ind(im,32); - % Write to the GIF File - if in == 1 - imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf); - else - imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY); - end + switch format + case '.gif' + im = frame2im(frame); + [imind,cm] = rgb2ind(im,32); + % Write to the GIF File + if in == 1 + imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf); + else + imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY); + end + case '.avi' + writeVideo(vidfile,frame); + otherwise + disp('Unknown format'); + break + end % terminal info while nbytes > 0 fprintf('\b') nbytes = nbytes - 1; end nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:))); in = in + 1; end disp(' ') - disp(['Gif saved @ : ',GIFNAME]) + switch format + case '.gif' + disp(['Gif saved @ : ',FILENAME]) + case '.avi' + disp(['Video saved @ : ',FILENAME]) + close(vidfile); + end +end end - diff --git a/matlab/create_mov.m b/matlab/create_mov.m deleted file mode 100644 index 5f84d65..0000000 --- a/matlab/create_mov.m +++ /dev/null @@ -1,69 +0,0 @@ -title1 = GIFNAME; -FIGDIR = BASIC.RESDIR; -GIFNAME = [FIGDIR, GIFNAME]; -XNAME = latexize(XNAME); -YNAME = latexize(YNAME); -FIELDNAME = latexize(FIELDNAME); - -vidfile = VideoWriter(GIFNAME,'Uncompressed AVI'); -vidfile.FrameRate = FPS; -open(vidfile); -% Set colormap boundaries -hmax = max(max(max(FIELD(:,:,FRAMES)))); -hmin = min(min(min(FIELD(:,:,FRAMES)))); -scale = -1; -flag = 0; -if hmax == hmin - disp('Warning : h = hmin = hmax = const') -else -% Setup figure frame -figure('Color','white','Position', [100, 100, 500, 500]); - pcolor(X,Y,FIELD(:,:,1)); % to set up -% colormap gray - colormap(bluewhitered) - axis tight manual % this ensures that getframe() returns a consistent size - if INTERP - shading interp; - end - in = 1; - nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); - for n = FRAMES % loop over selected frames - scale = max(max(abs(FIELD(:,:,n)))); - pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot - if INTERP - shading interp; - end - if NORMALIZED - pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\ - caxis([-1,1]); - else - pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\ - if CONST_CMAP - caxis([-1,1]*max(abs([hmin hmax]))); % const color map - else - caxis([-1,1]*scale); % adaptive color map - end - title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... - ,', scale = ',sprintf('%.1e',scale)]); - end - set(pclr, 'edgecolor','none'); axis square; - xlabel(XNAME); ylabel(YNAME); %colorbar; - title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... - ,', scaling = ',sprintf('%.1e',scale)]); - drawnow - % Capture the plot as an image - frame = getframe(gcf); - writeVideo(vidfile,frame); - % terminal info - while nbytes > 0 - fprintf('\b') - nbytes = nbytes - 1; - end - nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:))); - in = in + 1; - end - disp(' ') - disp(['Video saved @ : ',GIFNAME]) - close(vidfile); -end - diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m index faee436..2fb7815 100644 --- a/matlab/load_marconi.m +++ b/matlab/load_marconi.m @@ -1,30 +1,31 @@ function [ RESDIR ] = load_marconi( outfilename ) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here hostfolder = outfilename; hostfile = [hostfolder,'out*']; MISCDIR = ['/misc/HeLaZ_outputs/',outfilename(46:end)]; RESDIR = ['../',outfilename(46:end)]; miscfolder = [MISCDIR,'.']; system(['mkdir -p ',miscfolder]); disp(['mkdir -p ',miscfolder]); resultfolder = [RESDIR,'.']; % CMD = ['rsync -r ahoffman@login.marconi.cineca.it:',hostfolder,'/* ',MISCDIR]; % disp(CMD); % system(CMD); % CMD = ['rsync -r --exclude ''outputs_*.h5'' ahoffman@login.marconi.cineca.it:',hostfolder,'/* ',MISCDIR]; % disp(CMD); % system(CMD); % SCP the output file from marconi to misc folder of SPCPC - CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder]; +% CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder]; + CMD = ['rsync ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder]; disp(CMD); system(CMD); % Load the fort.90 as well in misc folder CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort*.90',' ',miscfolder]; disp(CMD); system(CMD); % Put it also in the result directory CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort*.90',' ',resultfolder]; disp(CMD); system(CMD); end diff --git a/matlab/load_params.m b/matlab/load_params.m index fe3e992..7042b6c 100644 --- a/matlab/load_params.m +++ b/matlab/load_params.m @@ -1,64 +1,63 @@ -CO = h5readatt(filename,'/data/input','CO'); -% K_N = h5readatt(filename,'/data/input','eta_n'); -% K_T = h5readatt(filename,'/data/input','eta_T'); -K_N = h5readatt(filename,'/data/input','K_n'); -K_T = h5readatt(filename,'/data/input','K_T'); -K_E = h5readatt(filename,'/data/input','K_E'); -Q0 = h5readatt(filename,'/data/input','q0'); -SHEAR = h5readatt(filename,'/data/input','shear'); -EPS = h5readatt(filename,'/data/input','eps'); -PMAXI = h5readatt(filename,'/data/input','pmaxi'); -JMAXI = h5readatt(filename,'/data/input','jmaxi'); -PMAXE = h5readatt(filename,'/data/input','pmaxe'); -JMAXE = h5readatt(filename,'/data/input','jmaxe'); -NON_LIN = h5readatt(filename,'/data/input','NON_LIN'); -NU = h5readatt(filename,'/data/input','nu'); -Nx = h5readatt(filename,'/data/input','Nx'); -Ny = h5readatt(filename,'/data/input','Ny'); -L = h5readatt(filename,'/data/input','Lx'); -CLOS = h5readatt(filename,'/data/input','CLOS'); -DT_SIM = h5readatt(filename,'/data/input','dt'); -MU = h5readatt(filename,'/data/input','mu'); +DATA.CO = h5readatt(filename,'/data/input','CO'); +DATA.K_N = h5readatt(filename,'/data/input','eta_n'); +DATA.K_T = h5readatt(filename,'/data/input','eta_T'); +% DATA.K_N = h5readatt(filename,'/data/input','K_n'); +% DATA.K_T = h5readatt(filename,'/data/input','K_T'); +% DATA.K_E = h5readatt(filename,'/data/input','K_E'); +DATA.Q0 = h5readatt(filename,'/data/input','q0'); +DATA.SHEAR = h5readatt(filename,'/data/input','shear'); +DATA.EPS = h5readatt(filename,'/data/input','eps'); +DATA.PMAXI = h5readatt(filename,'/data/input','pmaxi'); +DATA.JMAXI = h5readatt(filename,'/data/input','jmaxi'); +DATA.PMAXE = h5readatt(filename,'/data/input','pmaxe'); +DATA.JMAXE = h5readatt(filename,'/data/input','jmaxe'); +DATA.NON_LIN = h5readatt(filename,'/data/input','NON_LIN'); +DATA.NU = h5readatt(filename,'/data/input','nu'); +DATA.Nx = h5readatt(filename,'/data/input','Nx'); +DATA.Ny = h5readatt(filename,'/data/input','Ny'); +DATA.L = h5readatt(filename,'/data/input','Lx'); +DATA.CLOS = h5readatt(filename,'/data/input','CLOS'); +DATA.DT_SIM = h5readatt(filename,'/data/input','dt'); +DATA.MU = h5readatt(filename,'/data/input','mu'); % MU = str2num(filename(end-18:end-14)); %bad... -W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y'; -W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y'; -W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y'; -W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y'; -W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y'; +DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y'; +DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y'; +DATA.W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y'; +DATA.W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y'; +DATA.W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y'; -if NON_LIN == 'y' - NON_LIN = 1; +if DATA.NON_LIN == 'y' + DATA.NON_LIN = 1; else - NON_LIN = 0; + DATA.NON_LIN = 0; end -switch abs(CO) - case 0; CONAME = 'LB'; - case 1; CONAME = 'DG'; - case 2; CONAME = 'SG'; - case 3; CONAME = 'PA'; - case 4; CONAME = 'FC'; - otherwise; CONAME ='UK'; +switch abs(DATA.CO) + case 0; DATA.CONAME = 'LB'; + case 1; DATA.CONAME = 'DG'; + case 2; DATA.CONAME = 'SG'; + case 3; DATA.CONAME = 'PA'; + case 4; DATA.CONAME = 'FC'; + otherwise; DATA.CONAME ='UK'; end -if (CO <= 0); CONAME = [CONAME,'DK']; -else; CONAME = [CONAME,'GK']; +if (DATA.CO <= 0); DATA.CONAME = [DATA.CONAME,'DK']; +else; DATA.CONAME = [DATA.CONAME,'GK']; end -if (CLOS == 0); CLOSNAME = 'Trunc.'; -elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; -elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; +if (DATA.CLOS == 0); DATA.CLOSNAME = 'Trunc.'; +elseif(DATA.CLOS == 1); DATA.CLOSNAME = 'Clos. 1'; +elseif(DATA.CLOS == 2); DATA.CLOSNAME = 'Clos. 2'; end -if (PMAXE == PMAXI) && (JMAXE == JMAXI) - degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; +if (DATA.PMAXE == DATA.PMAXI) && (DATA.JMAXE == DATA.JMAXI) + degngrad = ['P_',num2str(DATA.PMAXE),'_J_',num2str(DATA.JMAXE)]; else - degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... - '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; + degngrad = ['Pe_',num2str(DATA.PMAXE),'_Je_',num2str(DATA.JMAXE),... + '_Pi_',num2str(DATA.PMAXI),'_Ji_',num2str(DATA.JMAXI)]; end degngrad = [degngrad,'_Kn_%1.1f_nu_%0.0e_',... - CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; -degngrad = sprintf(degngrad,[K_N,NU,MU]); -if ~NON_LIN; degngrad = ['lin_',degngrad]; end -resolution = [num2str(Nx),'x',num2str(Ny/2),'_']; -gridname = ['L_',num2str(L),'_']; -PARAMS = [resolution,gridname,degngrad]; -% BASIC.RESDIR = [SIMDIR,PARAMS,'/']; + DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e']; +degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MU]); +if ~DATA.NON_LIN; degngrad = ['lin_',degngrad]; end +resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_']; +gridname = ['L_',num2str(DATA.L),'_']; +DATA.PARAMS = [resolution,gridname,degngrad]; \ No newline at end of file diff --git a/matlab/photomaton.m b/matlab/photomaton.m new file mode 100644 index 0000000..f054b94 --- /dev/null +++ b/matlab/photomaton.m @@ -0,0 +1,24 @@ +function [ FIGURE ] = photomaton( DATA,OPTIONS ) +%UNTITLED5 Summary of this function goes here +%% Processing +toplot = process_field(DATA,OPTIONS); +FNAME = toplot.FILENAME; +FRAMES = toplot.FRAMES; + +% +TNAME = []; +FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 numel(FRAMES) 1]) + for i_ = 1:numel(FRAMES) + subplot(1,numel(FRAMES),i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))]; + pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))); set(pclr, 'edgecolor','none');pbaspect(toplot.ASPECT) + colormap(bluewhitered); + xlabel(toplot.XNAME); ylabel(toplot.YNAME);set(gca,'ytick',[]); + title(sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_)))); + if OPTIONS.INTERP + shading interp; + end + end + legend(['$',toplot.FIELDNAME,'$']); + FIGURE.FIGNAME = [FNAME,'_snaps',TNAME]; +end + diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m index a8cddde..597952e 100644 --- a/matlab/plots/plot_radial_transport_and_shear.m +++ b/matlab/plots/plot_radial_transport_and_shear.m @@ -1,73 +1,65 @@ +function [FIGURE] = plot_radial_transport_and_shear(DATA, TAVG_0, TAVG_1) %Compute steady radial transport tend = TAVG_1; tstart = TAVG_0; -[~,its0D] = min(abs(Ts0D-tstart)); -[~,ite0D] = min(abs(Ts0D-tend)); -SCALE = (2*pi/Nx/Ny)^2; -gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE; -gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE; -if W_HF - q_infty_avg = mean(HFLUX_X(its0D:ite0D))*SCALE; -end -% Compute steady shearing rate -tend = TAVG_1; tstart = TAVG_0; -[~,its2D] = min(abs(Ts3D-tstart)); -[~,ite2D] = min(abs(Ts3D-tend)); -shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1)); -shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1)); -% plots -fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) +[~,its0D] = min(abs(DATA.Ts0D-tstart)); +[~,ite0D] = min(abs(DATA.Ts0D-tend)); +SCALE = (2*pi/DATA.Nx/DATA.Ny)^2; +gamma_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE; +gamma_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE; + +FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600]) subplot(311) - yyaxis left - plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on; - plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',... +% yyaxis left + plot(DATA.Ts0D,DATA.PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on; + plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',... 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]); grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$') - ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]); - title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\kappa_N=$',num2str(K_N),... - ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... - ' $\mu_{hd}=$',num2str(MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]); - if W_HF - yyaxis right - plot(Ts0D,HFLUX_X*SCALE,'DisplayName','$\langle T_i \partial_y\phi \rangle_y$'); hold on; - grid on; set(gca,'xticklabel',[]); ylabel('$Q_x$') - xlim([Ts0D(1),Ts0D(end)]); ylim([0,5*abs(q_infty_avg)]); + ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]); + title(['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ', $\kappa_N=$',num2str(DATA.K_N),... + ', $L=',num2str(DATA.L),'$, $N=',num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',num2str(DATA.JMAXI),')$,',... + ' $\mu_{hd}=$',num2str(DATA.MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]); + %% zonal vs nonzonal energies for phi(t) and shear radial profile + + % computation + Ns3D = numel(DATA.Ts3D); + [KY, KX] = meshgrid(DATA.ky, DATA.kx); + Ephi_Z = zeros(1,Ns3D); + Ephi_NZ_kgt0 = zeros(1,Ns3D); + Ephi_NZ_kgt1 = zeros(1,Ns3D); + Ephi_NZ_kgt2 = zeros(1,Ns3D); + high_k_phi = zeros(1,Ns3D); + dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D); +% dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D); + for it = 1:numel(DATA.Ts3D) + [amp,ikzf] = max(abs((KX~=0).*DATA.PHI(:,1,1,it))); + Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2)))); + Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2)))); + Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2)))); + Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(DATA.PHI(:,:,1,it)).^2)))); + dxphi(:,:,1,it) = real(fftshift(ifft2(-1i*KX.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny))); +% dx2phi(:,:,1,it) = real(fftshift(ifft2(-KX.^2.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny))); end - % -% subplot(312) -% clr = line_colors(1,:); -% lstyle = line_styles(1); -% % plt = @(x_) mean(x_,1); -% plt = @(x_) x_(1,:); -% plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(\partial^2_x\phi)$'); hold on; -% plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle \partial^2_x\phi\rangle_y$'); hold on; -% plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle \partial^2_x\phi\rangle_x$'); hold on; -% plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle \partial^2_x\phi\rangle_{x,y}$'); hold on; -% plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',... -% 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]); -% ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]); -% grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show'); - %% zonal vs nonzonal energies for phi(t) + [~,its3D] = min(abs(DATA.Ts3D-tstart)); + [~,ite3D] = min(abs(DATA.Ts3D-tend)); + % plot subplot(312) it0 = 1; itend = Ns3D; trange = it0:itend; pltx = @(x) x;%-x(1); - plty = @(x) x./max((x(1:end))); + plty = @(x) x./max((x(its3D:ite3D))); % plty = @(x) x(500:end)./max(squeeze(x(500:end))); yyaxis left - plot(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on; - ylim([1e-1, 1.5]); ylabel('$E_{Z}$') + plot(pltx(DATA.Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',DATA.CONAME]); hold on; + ylim([-0.1, 1.5]); ylabel('$E_{Z}$') yyaxis right - plot(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$'); -% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName','$k_p>1$'); -% semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName','$k_p>2$'); - % semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]); -% legend('Location','southeast') - xlim([Ts3D(it0), Ts3D(itend)]); - ylim([1e-6, 1.0]); ylabel('$E_{NZ}$') + plot(pltx(DATA.Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$'); + xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]); + ylim([-0.1, 1.5]); ylabel('$E_{NZ}$') xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]); subplot(313) - [TY,TX] = meshgrid(x,Ts3D); - pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar; + [TY,TX] = meshgrid(DATA.x,DATA.Ts3D); + pclr = pcolor(TX,TY,squeeze((mean(dxphi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x\phi\rangle_y$') %colorbar; +% pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar; caxis([-3 3]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256)) -save_figure +end \ No newline at end of file diff --git a/matlab/post_processing.m b/matlab/post_processing.m index f30b11d..1d051ec 100644 --- a/matlab/post_processing.m +++ b/matlab/post_processing.m @@ -1,225 +1,205 @@ %% Retrieving max polynomial degree and sampling info Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe); Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi); Ns5D = numel(Ts5D); Ns3D = numel(Ts3D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts3D = Ts3D'; -%% Build grids -Nkx = numel(kx); Nky = numel(ky); -[KY,KX] = meshgrid(ky,kx); -Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky); -dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1); -KPERP2 = KY.^2+KX.^2; -[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky)); -[KY_XY,KX_XY] = meshgrid(ky,kx); -[KZ_XZ,KX_XZ] = meshgrid(z,kx); -[KZ_YZ,KY_YZ] = meshgrid(z,ky); - -Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z); -Lx = 2*pi/dkx; Ly = 2*pi/dky; -dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz; -x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x); -y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y); -[Y_XY,X_XY] = meshgrid(y,x); -[Z_XZ,X_XZ] = meshgrid(z,x); -[Z_YZ,Y_YZ] = meshgrid(z,y); - %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density ni00 = zeros(Nx,Ny,Nz,Ns3D); % " phi = zeros(Nx,Ny,Nz,Ns3D); % Electrostatic potential Z_phi = zeros(Nx,Ny,Nz,Ns3D); % Zonal " dens_e = zeros(Nx,Ny,Nz,Ns3D); % Particle density dens_i = zeros(Nx,Ny,Nz,Ns3D); %" Z_n_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal " Z_n_i = zeros(Nx,Ny,Nz,Ns3D); %" temp_e = zeros(Nx,Ny,Nz,Ns3D); % Temperature temp_i = zeros(Nx,Ny,Nz,Ns3D); % " Z_T_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal " Z_T_i = zeros(Nx,Ny,Nz,Ns3D); %" dyTe = zeros(Nx,Ny,Nz,Ns3D); % Various derivatives dyTi = zeros(Nx,Ny,Nz,Ns3D); % " dyni = zeros(Nx,Ny,Nz,Ns3D); % " dxphi = zeros(Nx,Ny,Nz,Ns3D); % " dyphi = zeros(Nx,Ny,Nz,Ns3D); % " dx2phi = zeros(Nx,Ny,Nz,Ns3D); % " for it = 1:numel(Ts3D) for iz = 1:numel(z) if KIN_E NE_ = Ne00(:,:,iz,it); ne00 (:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny))); if(W_DENS) DENS_E_ = DENS_E(:,:,iz,it); dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny))); Z_n_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_.*(KY==0)),Nx,Ny))); end if(W_TEMP) TEMP_E_ = TEMP_E(:,:,iz,it); dyTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny))); temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny))); Z_T_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny))); end end NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); ni00 (:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny))); phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny))); Z_phi (:,:,iz,it) = real(fftshift(ifft2((PH_.*(KY==0)),Nx,Ny))); dxphi (:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny))); dx2phi(:,:,iz,it) = real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny))); dyphi (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny))); if(W_DENS) DENS_I_ = DENS_I(:,:,iz,it); dyni (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny))); dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny))); Z_n_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_.*(KY==0)),Nx,Ny))); end if(W_TEMP) TEMP_I_ = TEMP_I(:,:,iz,it); dyTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny))); temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny))); Z_T_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_.*(KY==0)),Nx,Ny))); end end end % Post processing disp('- post processing') % particle flux Gamma_x= zeros(Nx,Ny,Nz,Ns3D); % Radial particle transport Gamma_y= zeros(Nx,Ny,Nz,Ns3D); % Azymuthal particle transport % heat flux Q_x = zeros(Nx,Ny,Nz,Ns3D); Q_y = zeros(Nx,Ny,Nz,Ns3D); % Epot averages phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % . % Kperp spectrum interpolation %full kperp points kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]); % interpolated kperps nk_noAA = floor(2/3*numel(kx)); kp_ip = kx; [thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip); [xn,yn] = pol2cart(thg,rg); [ky_s, sortIdx] = sort(ky); [xc,yc] = meshgrid(ky_s,kx); phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D); % for it = 1:numel(Ts3D) % Loop over 2D aX_XYays for iz = 1:numel(z) phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it)))); phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it)))); phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it)))); phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it)))); if(W_DENS) Gamma_x(:,:,iz,it) = -dens_i(:,:,iz,it).*dyphi(:,:,iz,it); Gamma_y(:,:,iz,it) = dens_i(:,:,iz,it).*dxphi(:,:,iz,it); end if(W_TEMP) Q_x(:,:,iz,it) = -temp_e(:,:,iz,it).*dyphi(:,:,iz,it); Q_y(:,:,iz,it) = temp_i(:,:,iz,it).*dxphi(:,:,iz,it); end shear_maxx_maxy(iz,it) = max( max(squeeze(-(dx2phi(:,:,iz,it))))); shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dx2phi(:,:,iz,it))))); shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dx2phi(:,:,iz,it))))); shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dx2phi(:,:,iz,it))))); Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn); phi_kp_t(:,iz,it) = mean(Z_rth,2); end end % for it = 1:numel(Ts5D) % Loop over 5D aX_XYays [~, it2D] = min(abs(Ts3D-Ts5D(it))); if KIN_E Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky; end Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky; end %% Compute primary instability growth rate % disp('- growth rate') % % Find max value of transport (end of linear mode) % [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2); % [~,itmax] = min(abs(Ts3D-tmax)); % tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax); % [~,its3D_lin] = min(abs(Ts3D-tstart)); % [~,ite3D_lin] = min(abs(Ts3D-tend)); % % g_I = zeros(Nkx,Nky,Nz); % for ikx = 1:Nkx % for iky = 1:Nky % for iz = 1:Nz % [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); % end % end % end % [gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3); % kmax_I = abs(ky(ikmax_I)); % Bohm_transport = K_N*gmax_I/kmax_I^2; %% Compute secondary instability growth rate % disp('- growth rate') % % Find max value of transport (end of linear mode) % % [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2); % % [~,itmax] = min(abs(Ts2D-tmax)); % % tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax); % [~,its3D_lin] = min(abs(Ts3D-tstart)); % [~,ite3D_lin] = min(abs(Ts3D-tend)); % % g_II = zeros(Nkx,Nky); % for ikx = 1:Nkx % for iky = 1 % for iz = 1:Nz % [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); % end % end % end % [gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3); % kmax_II = abs(kx(ikmax_II)); %% zonal vs nonzonal energies for phi(t) Ephi_Z = zeros(1,Ns3D); Ephi_NZ_kgt0 = zeros(1,Ns3D); Ephi_NZ_kgt1 = zeros(1,Ns3D); Ephi_NZ_kgt2 = zeros(1,Ns3D); high_k_phi = zeros(1,Ns3D); for it = 1:numel(Ts3D) % Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2))); % Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2))); [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it))); % Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))); Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))); Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))); Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))); % Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2; Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2)))); % Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it); % high_k_phi(it) = squeeze(abs(PHI(18,18,1,it)).^2); % kperp sqrt(2) % high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5 end diff --git a/matlab/process_field.m b/matlab/process_field.m new file mode 100644 index 0000000..563694e --- /dev/null +++ b/matlab/process_field.m @@ -0,0 +1,196 @@ +function [ TOPLOT ] = process_field( DATA, OPTIONS ) +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +%% Setup time +%% +FRAMES = zeros(size(OPTIONS.TIME)); +for i = 1:numel(OPTIONS.TIME) + [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D)); +end +%% Setup the plot geometry +[KY,~] = meshgrid(DATA.ky,DATA.kx); +directions = {'x','y','z'}; +Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(OPTIONS.TIME); +POLARPLOT = OPTIONS.POLARPLOT; +LTXNAME = OPTIONS.NAME; +switch OPTIONS.PLAN + case 'xy' + XNAME = '$x$'; YNAME = '$y$'; + [Y,X] = meshgrid(DATA.y,DATA.x); + REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; + case 'xz' + XNAME = '$x$'; YNAME = '$z$'; + [Y,X] = meshgrid(DATA.z,DATA.x); + REALP = 1; COMPDIM = 2; SCALE = 0; + case 'yz' + XNAME = '$y$'; YNAME = '$z$'; + [Y,X] = meshgrid(DATA.z,DATA.y); + REALP = 1; COMPDIM = 1; SCALE = 0; + case 'kxky' + XNAME = '$k_x$'; YNAME = '$k_y$'; + [Y,X] = meshgrid(DATA.ky,DATA.kx); + REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; + case 'kxz' + XNAME = '$k_x$'; YNAME = '$z$'; + [Y,X] = meshgrid(DATA.z,DATA.kx); + REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0; + OPTIONS.COMP = 'max'; + case 'kyz' + XNAME = '$k_y$'; YNAME = '$z$'; + [Y,X] = meshgrid(DATA.z,DATA.ky); + REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0; + OPTIONS.COMP = 'max'; +end +Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y))); +Lmin_ = min([Xmax_,Ymax_]); +Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2; +Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2; +DIMENSIONS = [100, 100, 400*Rx, 400*Ry]; +ASPECT = [Rx, Ry, 1]; +% Polar grid +POLARNAME = []; +if POLARPLOT + POLARNAME = 'polar'; + X__ = (X+DATA.a).*cos(Y); + Y__ = (X+DATA.a).*sin(Y); + X = X__; + Y = Y__; + XNAME='X'; + YNAME='Z'; + DIMENSIONS = [100, 100, 500, 500]; + ASPECT = [1,1,1]; + sz = size(X); + FIELD = zeros(sz(1),sz(2),Nt); +else + sz = size(X); + FIELD = zeros(sz(1),sz(2),Nt); +end +%% Process the field to plot + +switch OPTIONS.COMP + case 'avg' + compr = @(x) mean(x,COMPDIM); + COMPNAME = ['avg',directions{COMPDIM}]; + FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}]; + case 'max' + compr = @(x) max(x,[],COMPDIM); + COMPNAME = ['max',directions{COMPDIM}]; + FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME]; + otherwise + switch COMPDIM + case 1 + i = OPTIONS.COMP; + compr = @(x) x(i,:,:); + COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i)); + FIELDNAME = [LTXNAME,'(',COMPNAME,')']; + case 2 + i = OPTIONS.COMP; + compr = @(x) x(:,i,:); + COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.y(i)); + FIELDNAME = [LTXNAME,'(',COMPNAME,')']; + case 3 + i = OPTIONS.COMP; + compr = @(x) x(:,:,i); + COMPNAME = sprintf([directions{COMPDIM},'=','%2.1f'],DATA.z(i)/pi); + FIELDNAME = [LTXNAME,'(',COMPNAME,')']; + end +end + +switch REALP + case 1 % Real space plot + process = @(x) real(fftshift(ifft2(x,Nx,Ny))); + shift_x = @(x) x; + shift_y = @(x) x; + case 0 % Frequencies plot + OPTIONS.INTERP = 0; + switch COMPDIM + case 1 + process = @(x) abs(fftshift(x,2)); + shift_x = @(x) fftshift(x,1); + shift_y = @(x) fftshift(x,1); + case 2 + process = @(x) abs((x)); + shift_x = @(x) (x); + shift_y = @(x) (x); + case 3 + process = @(x) abs(fftshift(x,2)); + shift_x = @(x) fftshift(x,2); + shift_y = @(x) fftshift(x,2); + end +end + +switch OPTIONS.NAME + case '\phi' + NAME = 'phi'; + if COMPDIM == 3 + for it = FRAMES + tmp = squeeze(compr(DATA.PHI(:,:,:,it))); + FIELD(:,:,it) = squeeze(process(tmp)); + end + else + if REALP + tmp = zeros(Nx,Ny,Nz); + else + tmp = zeros(DATA.Nkx,DATA.Nky,Nz); + end + for it = FRAMES + for iz = 1:numel(DATA.z) + tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,it))); + end + FIELD(:,:,it) = squeeze(compr(tmp)); + end + end + case 'n_i' + NAME = 'ni'; + if COMPDIM == 3 + for it = FRAMES + tmp = squeeze(compr(DATA.DENS_I(:,:,:,it))); + FIELD(:,:,it) = squeeze(process(tmp)); + end + else + if REALP + tmp = zeros(Nx,Ny,Nz); + else + tmp = zeros(DATA.Nkx,DATA.Nky,Nz); + end + for it = FRAMES + for iz = 1:numel(DATA.z) + tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it))); + end + FIELD(:,:,it) = squeeze(compr(tmp)); + end + end + case 'n_i^{NZ}' + NAME = 'niNZ'; + if COMPDIM == 3 + for it = FRAMES + tmp = squeeze(compr(DATA.DENS_I(:,:,:,it).*(KY~=0))); + FIELD(:,:,it) = squeeze(process(tmp)); + end + else + if REALP + tmp = zeros(Nx,Ny,Nz); + else + tmp = zeros(DATA.Nkx,DATA.Nky,Nz); + end + for it = FRAMES + for iz = 1:numel(DATA.z) + tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,it).*(KY~=0))); + end + FIELD(:,:,it) = squeeze(compr(tmp)); + end + end +end +TOPLOT.FIELD = FIELD; +TOPLOT.X = shift_x(X); +TOPLOT.Y = shift_y(Y); +TOPLOT.FIELDNAME = FIELDNAME; +TOPLOT.XNAME = XNAME; +TOPLOT.YNAME = YNAME; +TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME]; +TOPLOT.DIMENSIONS= DIMENSIONS; +TOPLOT.ASPECT = ASPECT; +TOPLOT.FRAMES = FRAMES; + +end + diff --git a/matlab/save_figure.m b/matlab/save_figure.m index c5a9c72..1aa1dd6 100644 --- a/matlab/save_figure.m +++ b/matlab/save_figure.m @@ -1,8 +1,10 @@ %% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME) % and parameters -if ~exist([BASIC.RESDIR,'/fig'], 'dir') - mkdir([BASIC.RESDIR,'/fig']) +function save_figure(data,FIGURE) +if ~exist([data.localdir,'/fig'], 'dir') + mkdir([data.localdir,'/fig']) end -saveas(fig,[BASIC.RESDIR,'/fig/', FIGNAME,'.fig']); -saveas(fig,[BASIC.RESDIR, FIGNAME,'.png']); -disp(['Figure saved @ : ',[BASIC.RESDIR, FIGNAME,'.png']]) \ No newline at end of file +saveas(FIGURE.fig,[data.localdir,'/fig/', FIGURE.FIGNAME,'.fig']); +saveas(FIGURE.fig,[data.localdir, FIGURE.FIGNAME,'.png']); +disp(['Figure saved @ : ',[data.localdir, FIGURE.FIGNAME,'.png']]) +end \ No newline at end of file diff --git a/matlab/show_geometry.m b/matlab/show_geometry.m new file mode 100644 index 0000000..e57746d --- /dev/null +++ b/matlab/show_geometry.m @@ -0,0 +1,96 @@ +function [ FIGURE ] = show_geometry(DATA,OPTIONS) +% filtering Z pinch and torus +if DATA.Nz > 1 % Torus flux tube geometry + Nptor = DATA.Nz; Tgeom = 1; + Nturns = 1; + a = DATA.EPS; % Torus minor radius + R = 1.; % Torus major radius + q = (DATA.Nz>1)*DATA.Q0; % Flux tube safety factor + theta = linspace(-pi, pi, Nptor) ; % Poloidal angle + phi = linspace(0., 2.*pi, Nptor) ; % Toroidal angle + [t, p] = meshgrid(phi, theta); + x_tor = (R + a.*cos(p)) .* cos(t); + y_tor = (R + a.*cos(p)) .* sin(t); + z_tor = a.*sin(p); +else % Zpinch geometry + Nturns = 0.1; Tgeom = 0; + a = 1; % Torus minor radius + R = 0; % Torus major radius + q = 0; + theta = linspace(-Nturns*pi, Nturns*pi, 100); % Poloidal angle + phi = linspace(-0.1, 0.1, 3); % cylinder height + [t, p] = meshgrid(phi, theta); + x_tor = a.*cos(p); + z_tor = a.*sin(p); + y_tor = t; +end +FIGURE.fig = figure; set(gcf, 'Position', [50 50 1200 600]) +torus=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local') +set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none') +xlabel('X');ylabel('Y');zlabel('Z'); +% field line coordinates +Xfl = @(z) (R+a*cos(z)).*cos(q*z); +Yfl = @(z) (R+a*cos(z)).*sin(q*z); +Zfl = @(z) a*sin(z); +Rvec= @(z) [Xfl(z); Yfl(z); Zfl(z)]; +% xvec +xX = @(z) (Xfl(z)-R*cos(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2); +xY = @(z) (Yfl(z)-R*sin(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2); +xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2); +xvec= @(z) [xX(z); xY(z); xZ(z)]; +% bvec +bX = @(z) Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z))./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2); +bY = @(z) (a*cos(z).*sin(q*z) + q*Xfl(z))./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2); +bZ = @(z) a*cos(z)./sqrt(Tgeom*(a*cos(z).*cos(q*z) - q*Yfl(z)).^2+(a*cos(z).*sin(q*z) + q*Xfl(z)).^2+(a*cos(z)).^2); +bvec= @(z) [bX(z); bY(z); bZ(z)]; +% yvec = b times x +yX = @(z) bY(z).*xZ(z)-bZ(z).*xY(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2); +yY = @(z) bZ(z).*xX(z)-bX(z).*xZ(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2); +yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2); +yvec= @(z) [yX(z); yY(z); yZ(z)]; +% Plot high res field line +theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle +plot3(Xfl(theta),Yfl(theta),Zfl(theta)) +% Planes plot +theta = DATA.z ; % Poloidal angle +plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok') +% Plot x,y,bvec at these points +scale = 0.15; +quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r'); +quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g'); +quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b'); + +%% Plot plane result +OPTIONS.POLARPLOT = 0; +OPTIONS.PLAN = 'xy'; +r_o_R = DATA.rho_o_R; +[Y,X] = meshgrid(r_o_R*DATA.y,r_o_R*DATA.x); +max_ = 0; +FIELDS = zeros(DATA.Nx,DATA.Ny,numel(OPTIONS.PLANES)); + +for it_ = 1:numel(OPTIONS.TIME) +[~,it] = min(abs(OPTIONS.TIME(it_)-DATA.Ts3D)); +for iz = OPTIONS.PLANES + OPTIONS.COMP = iz; + toplot = process_field(DATA,OPTIONS); + FIELDS(:,:,iz) = toplot.FIELD(:,:,it); + tmp = max(max(abs(FIELDS(:,:,iz)))); + if (tmp > max_) max_ = tmp; +end +for iz = OPTIONS.PLANES + z_ = DATA.z(iz); + Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_); + Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_); + Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_); + s=surface(Xp,Yp,Zp,FIELDS(:,:,iz)/max_); + s.EdgeColor = 'none'; + colormap(bluewhitered); + caxis([-1,1]); +end +end +%% +axis equal +view([1,-2,1]) + +end + diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index 5515dca..249c2a4 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,353 +1,213 @@ addpath(genpath('../matlab')) % ... add addpath(genpath('../matlab/plots')) % ... add outfile =''; %% Directory of the simulation -if 1% Local results +if 0% Local results outfile =''; -outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_adiabe'; -% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_lin'; -% outfile ='fluxtube_salphaB_s0/64x64x16_5x3_L_200_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK'; -% outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK'; -% outfile ='Linear_Device/64x64x20_5x2_Lx_20_Ly_150_q0_25_kN_0.24_kT_0.03_nu_1e-02_DGGK'; - BASIC.RESDIR = ['../results/',outfile,'/']; - BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; - system(['mkdir -p ',BASIC.MISCDIR]); - CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD); +outfile =''; +outfile =''; +outfile =''; +outfile ='fluxtube_salphaB_s0/100x100x30_5x3_Lx_200_Ly_100_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe'; +% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe'; +% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_8_nu_1e-01_DGGK_adiabe_Sg'; + LOCALDIR = ['../results/',outfile,'/']; + MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; + system(['mkdir -p ',MISCDIR]); + CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); else% Marconi results outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x200_L_200_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt'; % outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt'; -BASIC.RESDIR = ['../',outfile(46:end-8),'/']; -BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/']; +% BASIC.RESDIR = ['../',outfile(46:end-8),'/']; +MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/']; end %% Load the results % Load outputs from jobnummin up to jobnummax -JOBNUMMIN = 00; JOBNUMMAX = 00; -% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A -compile_results %Compile the results from first output found to JOBNUMMAX if existing +JOBNUMMIN = 00; JOBNUMMAX = 20; +data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing -%% Post-processing -post_processing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) -TAVG_0 = 1500; TAVG_1 = 4000; % Averaging times duration -plot_radial_transport_and_shear +TAVG_0 = 1500; TAVG_1 = 5000; % Averaging times duration +fig = plot_radial_transport_and_shear(data,TAVG_0,TAVG_1); +save_figure(data,fig) end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options -t0 =000; iz = 1; ix = 1; iy = 1; -skip_ =1; FPS = 30; DELAY = 1/FPS; -[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D); -[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D); -T = Ts3D; FRAMES = FRAMES_3D; -INTERP = 0; NORMALIZED = 1; CONST_CMAP = 0; BWR =1;% Gif options -% Field to plot -% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i'; -% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}'; -% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i'; -% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}'; -% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}'; -% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}'; -FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi'; -% FIELD = phi-Z_phi; NAME = 'NZphi'; FIELDNAME = '\phi_{NZ}'; -% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x'; - -% Sliced -% plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; -% plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; -plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; - -% K-space -% FIELD = PHI.*(KY~=0); NAME = 'NZPHI'; FIELDNAME = '\tilde \phi_{k_y\neq0}'; -% FIELD = Ne00; NAME = 'Ne00'; FIELDNAME = 'N_e^{00}'; -% FIELD = Ni00; NAME = 'Ni00'; FIELDNAME = 'N_i^{00}'; -% plt = @(x) fftshift((abs(x( :, :,1,:))),2); X = fftshift(KX,2); Y = fftshift(KY,2); XNAME = 'k_x'; YNAME = 'k_y'; - -% Averaged -% plt = @(x) mean(x,1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; -% plt = @(x) mean(x,2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; -% plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; - -FIELD = squeeze(plt(FIELD)); - -% Naming -GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS]; - -% Create movie (gif or mp4) -create_gif -% create_mov +options.INTERP = 1; +options.POLARPLOT = 0; +options.NAME = 'n_i^{NZ}'; +options.PLAN = 'xy'; +options.COMP = 1; +options.TIME = 00:20:6000; +% options.TIME = 140:0.5:160; +data.a = data.EPS * 2000; +create_film(data,options,'.gif') end if 0 -%% 2D plot : real space - -% Chose the field to plot -% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}'; -% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}' -% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i'; -% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e'; -% FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}'; -% FIELD = dens_i-Z_n_i; FNAME = 'ni_NZ'; FIELDLTX = 'n_i^{NZ}'; -% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i'; -% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e'; -% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi'; -FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}'; -% FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}'; -% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x'; -% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}'; - -% Chose when to plot it -tf = [0 15 27 28 30]; - -% Planar plot: choose a plane to plot at x0/y0/z0 coordinates -x0 = 0.0; y0 = 0.0; z0 = 0.0; -[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0)); -% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')'] -% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')'] -plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(z=',num2str(round(z(iz)/pi)),'\pi)'] - -% Averaged -% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_x'] -% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_y'] -% plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_z'] - - -% -TNAME = []; -fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350]) -plt_2 = @(x) x;%./max(max(x)); - for i_ = 1:numel(tf) - [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; - subplot(1,numel(tf),i_) - DATA = plt_2(squeeze(plt(FIELD,it))); - pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) - colormap(bluewhitered); %caxis([-30,30]); - xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]); - title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); - end - legend(latexize(FIELDLTX)); -save_figure +%% 2D snapshots +% Options +options.INTERP = 0; +options.POLARPLOT = 0; +options.NAME = 'n_i^{NZ}'; +options.PLAN = 'xy'; +options.COMP = 1; +options.TIME = [100 400 2500 5500]; +data.a = data.EPS * 1000; +fig = photomaton(data,options); +save_figure(data,fig) end if 0 -%% 2D plot : k space - -% Chose the field to plot -% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}'; -FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}' -% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi'; -% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x'; -% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e'; - -% Chose when to plot it -% tf = [0 15 27 28 30]; -tf = 200:200:1200; -% tf = 8000; - -% Planar plot: choose a plane to plot at x0/y0/z0 coordinates -x0 = 0.0; y0 = 0.3; z0 = 0.5*pi; -[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0)); -% plt = @(x,it) abs(x(ix, :, :,it)); X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_x=',num2str(round(kx(ix))),')']; -% plt = @(x,it) abs(x( :,iy, :,it)); X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_y=',num2str(round(ky(iy))),')']; -% plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)']; -% % -% plt_x = @(x) fftshift(x,1); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,1); plt = @(x,it) max(abs(x( :, :,:,it)),[],1); -% X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_x)']; -% -% plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],2); -% X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_y)']; - -plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,2); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],3); -X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(\max_z)']; - -% -TNAME = []; -fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 300*numel(tf), 500]) - for i_ = 1:numel(tf) - [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; - subplot(1,numel(tf),i_) - DATA = plt_2(squeeze(plt(FIELD,it))); - pclr = pcolor(plt_x(X),plt_y(Y),plt_z(DATA)); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1]) - caxis([0 1]*1e3); -% caxis([-1 1]*5); - xlabel(latexize(XNAME)); ylabel(latexize(YNAME)); - if(i_ > 1); set(gca,'ytick',[]); end; - title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); - end - legend(latexize(FIELDLTX)); -save_figure +%% 3D plot on the geometry +options.INTERP = 1; +options.NAME = 'n_i'; +options.PLANES = 1; +options.TIME = 5000; +data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio +FIGURE = show_geometry(data,options); end if 0 %% TAVG_0 = 1000; TAVG_1 = 5000; % Averaging times duration ZF_fourier_analysis end if 0 %% plot_param_evol end if 0 %% Radial shear profile (with moving average) tf = 1000+[0:100:1000]; ymax = 0; figure for i_ = 1:numel(tf) [~,it] = min(abs(Ts3D-tf(i_))); data = squeeze((mean(dx2phi(:,:,1,it),2))); step = 50; plot(movmean(x,step),movmean(data,step),'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on; ymax = max([ymax abs(min(data)) abs(max(data))]); end xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax); xlabel('$x/\rho_s$'); ylabel('$s_{E\times B,x}$'); grid on end if 0 %% zonal vs nonzonal energies for phi(t) t0 = 200; [~, it0] = min(abs(Ts3D-t0)); itend = Ns3D; trange = it0:itend; pltx = @(x) x;%-x(1); plty = @(x) x./max(squeeze(x)); fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS]; set(gcf, 'Position', [100, 100, 1400, 500]) subplot(121) % yyaxis left semilogy(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on; % yyaxis right semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName',['NZ, $k_p>0$, ',CONAME]); semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName',['NZ, $k_p>1$, ',CONAME]); semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName',['NZ, $k_p>2$, ',CONAME]); % semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]); title('Energy'); legend('Location','southeast') xlim([Ts3D(it0), Ts3D(itend)]); ylim([1e-3, 1.5]) xlabel('$t c_s/R$'); grid on;% xlim([0 500]); subplot(122) % plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange))); plot3(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)),Ts3D(trange)); title('Phase space'); legend(CONAME) xlabel('$E_Z$'); ylabel('$E_{NZ}$'); zlabel('time'); grid on;% xlim([0 500]); end if 0 %% Conservation laws Nxmax = Nx/2; Nymax = Ny/2; mflux_x_i = squeeze(sum((Gamma_x( 1,1:Nxmax,:)+Gamma_x( 1,2:Nxmax+1,:))/2,2)./sum(Gamma_x( 1,1:Nxmax,:))); mflux_x_o = squeeze(sum((Gamma_x( Nxmax,1:Nxmax,:)+Gamma_x( Nxmax,2:Nxmax+1,:))/2,2)./sum(Gamma_x( Nxmax,1:Nxmax,:))); mflux_y_i = squeeze(sum((Gamma_y(1:Nxmax, 1,:)+Gamma_y(2:Nxmax+1, 1,:))/2,1)./sum(Gamma_y(1:Nxmax, 1,:))); mflux_y_o = squeeze(sum((Gamma_y(1:Nxmax, Nymax,:)+Gamma_y(2:Nxmax+1, Nymax,:))/2,1)./sum(Gamma_y(1:Nxmax, Nymax,:))); mass_cons = mflux_x_i - mflux_x_o + mflux_y_i - mflux_y_o; %% figure plt = @(x) squeeze(mean(mean(x(:,:,1,:),1),2)); subplot(211) plot(Ts3D,plt(dens_e+dens_i),'DisplayName','$\delta n_e + \delta n_i$'); hold on; plot(Ts3D,plt(ne00+ni00),'DisplayName','$\delta n_e^{00} + \delta n_i^{00}$'); hold on; plot(Ts3D,plt(temp_e+temp_i),'DisplayName','$\delta T_e + \delta T_i$'); hold on; legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]) subplot(212); plot(Ts3D,mass_cons*(2*pi/Nx/Ny)^2,'DisplayName','in-out'); hold on % plot(Ts3D,squeeze(mflux_x_i),'DisplayName','Flux i x'); % plot(Ts3D,squeeze(mflux_x_o),'DisplayName','Flux o x'); % plot(Ts3D,squeeze(mflux_y_i),'DisplayName','Flux i y'); % plot(Ts3D,squeeze(mflux_y_o),'DisplayName','Flux o y'); legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]); %ylim([-0.1, 2]*mean(mflux_x_i)) end if 0 %% Zonal profiles (ky=0) % Chose the field to plot FIELD = Ne00.*conj(Ne00); FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2' % FIELD = Ni00.*conj(Ni00); FNAME = 'Ni002'; FIELDLTX = '|N_i^{00}|^2' % FIELD = abs(PHI); FNAME = 'absPHI'; FIELDLTX = '|\tilde\phi|'; % FIELD = PHI.*conj(PHI); FNAME = 'PHI2'; FIELDLTX = '|\tilde\phi|^2'; % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x'; % FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e'; % Chose when to plot it tf = 200:200:1200; % tf = 8000; % Sliced plt = @(x,it) x( 2:end, 1,1,it)./max(max(x( 2:end, 1,1,:))); XNAME = 'k_x'; % TNAME = []; fig = figure; FIGNAME = ['Zonal_',FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 600,400]) plt_2 = @(x) (fftshift(x,2)); for i_ = 1:numel(tf) [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; DATA = plt_2(squeeze(plt(FIELD,it))); semilogy(kx(2:end),DATA,'-','DisplayName',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on; grid on; xlabel(latexize(XNAME)); end title(['$',FIELDLTX,'$ Zonal Spectrum']); legend('show'); save_figure end if 0 %% Time evolutions and growth rate plot_time_evolution_and_gr end if 0 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3) % tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window tstart = 0000; tend = 1000; % tstart = 10000; tend = 12000; % Chose the field to plot % FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}'; % FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}' FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi'; % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x'; LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0; plot_kperp_spectrum -end - -if 0 -%% Torus plot -aminor = EPS; % Torus minor radius -Rmajor = 1.; % Torus major radius -theta = linspace(-pi, pi, 30) ; % Poloidal angle -phi = linspace(0., 2.*pi, 30) ; % Toroidal angle -[t, p] = meshgrid(phi, theta); -x = (Rmajor + aminor.*cos(p)) .* cos(t); -y = (Rmajor + aminor.*cos(p)) .* sin(t); -z = aminor.*sin(p); -figure; -torus=surf(x, y, z); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local') -set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none') -xlabel('X');ylabel('Y');zlabel('Z'); -% field line plot -Nturns = 1; -theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle -xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0); -yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0); -zFL = aminor.*sin(theta); -plot3(xFL,yFL,zFL,'r') -% Planes plot -theta = linspace(-pi, pi, Nz) ; % Poloidal angle -xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0); -yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0); -zFL = aminor.*sin(theta); -plot3(xFL,yFL,zFL,'sb') -axis equal - end \ No newline at end of file diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m index 95365b4..7f0d2c3 100644 --- a/wk/linear_1D_entropy_mode.m +++ b/wk/linear_1D_entropy_mode.m @@ -1,155 +1,156 @@ RUN = 1; % To run or just to load 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 TAU = 1.0; % e/i temperature ratio K_N = 1/0.6; % Density gradient drive K_T = 0.0; % Temperature ''' K_E = 0.0; % Electrostat ''' SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) %% GRID PARAMETERS -NX = 100; % real space x-gridpoints +NX = 150; % real space x-gridpoints NY = 1; % '' y-gridpoints -LX = 150; % Size of the squared frequency domain +LX = 200; % Size of the squared frequency domain LY = 1; % Size of the squared frequency domain NZ = 1; % number of perpendicular planes (parallel grid) Q0 = 1.0; % safety factor SHEAR = 0.0; % magnetic shear EPS = 0.0; % inverse aspect ratio +SG = 1; % Staggered z grids option %% TIME PARMETERS -TMAX = 100; % Maximal time unit +TMAX = 200; % Maximal time unit DT = 1e-2; % 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; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS SIMID = 'Linear_entropy_mode'; % Name of the simulation NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1) KIN_E = 1; % Collision operator % (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK) -CO = 2; +CO = 4; INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax)) NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_PHI= 0; % Start simulation with a noisy phi %% 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.5; % Hyper diffusivity cutoff ratio kmax = NX*pi/LX;% Highest fourier mode NU_HYP = 0.0; % Hyperdiffusivity coefficient MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 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 LAMBDAD = 0.0; NOISE0 = 1.0e-5; % Init noise amplitude BCKGD0 = 0.0; % Init background GRADB = 1.0; CURVB = 1.0; %% PARAMETER SCANS if 1 %% Parameter scan over PJ % PA = [2 4]; % JA = [1 2]; PA = [4]; JA = [2]; DTA= DT*ones(size(JA));%./sqrt(JA); % DTA= DT; mup_ = MU_P; muj_ = MU_J; Nparam = numel(PA); param_name = 'PJ'; gamma_Ni00 = zeros(Nparam,floor(NX/2)+1); gamma_Nipj = zeros(Nparam,floor(NX/2)+1); gamma_phi = zeros(Nparam,floor(NX/2)+1); for i = 1:Nparam % Change scan parameter PMAXE = PA(i); PMAXI = PA(i); JMAXE = JA(i); JMAXI = JA(i); DT = DTA(i); setup system(['rm fort*.90']); % Run linear simulation if RUN -% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk']) + system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk']) - system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk']) +% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk']) end % Load and process results %% filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results for ikx = 1:NX/2+1 tend = max(Ts3D(abs(Ni00(ikx,1,1,:))~=0)); tstart = 0.6*tend; [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); trange = itstart:itend; % exp fit on moment 00 X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,1,1,trange))); gamma_Ni00(i,ikx) = LinearFit_s(X_,Y_); % exp fit on phi X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,1,1,trange))); gamma_phi (i,ikx) = LinearFit_s(X_,Y_); end gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0)); gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0)); if 0 %% Fit verification figure; for i = 1:1:Nx/2+1 X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:))); semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on; end 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*kx),plt(gamma_phi(i,1:end)),... 'Color',clr,... 'LineStyle',linestyle{1},'Marker','^',... ...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']); 'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']); hold on; end grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(kx)]); title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$']) legend('show'); %xlim([0.01,10]) saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']); saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']); end end if 0 %% Space time [YT,XT] = meshgrid(Ts3D,kx); figure; % pclr = surf(XT,YT,squeeze(abs(PHI_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar; % pclr = pcolor(XT,YT,squeeze(abs(Ni00_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar; semilogy(Ts3D(1:TMAX/SPS3D),squeeze(abs(PHI_ST(1,50:5:100,:)))); end end \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index 10fefe6..5490687 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,81 +1,81 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency K_N = 2.22; % Density gradient drive K_T = 8.0; % Temperature ''' K_E = 0.00; % Electrostat gradient SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) NU_HYP = 0.0; KIN_E = 0; % Kinetic (1) or adiabatic (2) electron model %% GRID PARAMETERS -NX = 50; % Spatial radial resolution ( = 2x radial modes) -LX = 300; % Radial window size +NX = 100; % Spatial radial resolution ( = 2x radial modes) +LX = 200; % Radial window size NY = 100; % Spatial azimuthal resolution (= azim modes) -LY = 300; % Azimuthal window size -NZ = 20; % number of perpendicular planes (parallel grid) +LY = 100; % Azimuthal window size +NZ = 30; % number of perpendicular planes (parallel grid) P = 4; J = 2; %% GEOMETRY PARAMETERS Q0 = 2.7; % safety factor SHEAR = 0.0; % magnetic shear EPS = 0.18; % inverse aspect ratio GRADB = 1.0; % Magnetic gradient CURVB = 1.0; % Magnetic curvature SG = 1; % Staggered z grids option %% TIME PARAMETERS -TMAX = 50; % Maximal time unit -DT = 2e-2; % Time step +TMAX = 200; % Maximal time unit +DT = 5e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1; % Sampling per time unit for 2D arrays SPS3D = 5; % Sampling per time unit for 3D arrays SPS5D = 1/200; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints/10 JOB2LOAD= -1; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK) CO = 1; CLOS = 0; % Closure model (0: =0 truncation) -NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) +NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) SIMID = 'fluxtube_salphaB_s0'; % Name of the simulation % SIMID = 'simulation_A'; % Name of the simulation % SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1) % INIT options INIT_PHI= 0; % Start simulation with a noisy phi (0= noisy moments 00) INIT_ZF = 0; ZF_AMP = 0.0; INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; %% 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 PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; kmax = NX*pi/LX;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; BCKGD0 = 0.0; % Init background TAU = 1.0; % e/i temperature ratio MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% Setup and file management setup system('rm fort*.90'); outfile = [BASIC.RESDIR,'out.txt']; disp(outfile); diff --git a/wk/marconi_run.m b/wk/marconi_run.m index 15674da..e6ea6a5 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,118 +1,122 @@ clear all; addpath(genpath('../matlab')) % ... add SUBMIT = 1; % To submit the job automatically CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain) % EXECNAME = 'helaz_dbg'; - EXECNAME = 'helaz_3.9'; -for K_N = [1/0.6] + EXECNAME = 'helaz_3.5'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS CLUSTER.PART = 'prod'; % dbg or prod % CLUSTER.PART = 'dbg'; CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end; CLUSTER.MEM = '128GB'; % Memory CLUSTER.JNAME = 'HeLaZ';% Job name NP_P = 2; % MPI processes along p NP_KX = 24; % MPI processes along kx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency -K_N = 1.0/0.6; % Density gradient drive (R/Ln) +K_N = 1.0/0.6; % Density gradient drive (R/Ln) NU_HYP = 0.0; +SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) %% GRID PARAMETERS -Nx = 300; % Realspace x-gridpoints -Ny = 300; % Realspace y-gridpoints -Lx = 120; % Size of the squared frequency domain -Ly = 120; % Size of the squared frequency domain -Nz = 1; % number of perpendicular planes (parallel grid) -q0 = 1.0; % q factor () -shear = 0.0; % magnetic shear -eps = 0.0; % inverse aspect ratio -P = 8; -J = 4; +NX = 300; % Realspace x-gridpoints +NY = 300; % Realspace y-gridpoints +LX = 120; % Size of the squared frequency domain +LY = 120; % Size of the squared frequency domain +NZ = 1; % number of perpendicular planes (parallel grid) +Q0 = 1.0; % q factor () +SHEAR = 0.0; % magnetic shear +EPS = 0.0; % inverse aspect ratio +P = 4; +J = 2; %% TIME PARAMETERS TMAX = 10000; % Maximal time unit -DT = 8e-3; % Time step +DT = 2e-3; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1; % Sampling per time unit for 2D arrays SPS3D = 1/2; % Sampling per time unit for 3D arrays -SPS5D = 1/100; % Sampling per time unit for 5D arrays -JOB2LOAD= 2; % start from t=0 if <0, else restart from outputs_$job2load +SPS5D = 1/50; % Sampling per time unit for 5D arrays +JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK) -CO = 3; +CO = 2; CLOS = 0; % Closure model (0: =0 truncation) -NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) +NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) % SIMID = 'test_chained_job'; % Name of the simulation SIMID = 'simulation_A'; % Name of the simulation % SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1) % INIT options INIT_ZF = 0; ZF_AMP = 0.0; INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; %% 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 +KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model +GRADB = 1.0; % Magnetic gradient +CURVB = 1.0; % Magnetic curvature +SG = 0; % Staggered z grids option PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; -kmax = Nx*pi/Lx;% Highest fourier mode +kmax = NX*pi/LX;% 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; BCKGD0 = 0.0; % Init background TAU = 1.0; % e/i temperature ratio -K_T = 0.0; % Temperature gradient +K_T = 0.0; % Temperature gradient +K_E = 0.0; % ES ''' INIT_PHI= 1; % Start simulation with a noisy phi and moments MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f % Compute processes distribution Ntot = NP_P * NP_KX; Nnodes = ceil(Ntot/48); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kx CLUSTER.CPUPT = '1'; % CPU per task %% Run file management scripts setup SBATCH_CMD = 'sbatch batch_script.sh\n'; write_sbash_marconi if(mod(NP_P*NP_KX,48)~= 0) disp('WARNING : unused cores (ntot cores must be a 48 multiple)'); end if(SUBMIT) [~,job_info_] = system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); disp(job_info_); jobid_ = job_info_(21:27); if(CHAIN>0) for CHAIN_IDX = 1:CHAIN SBATCH_CMD = ['sbatch --dependency=afterok:',jobid_,' batch_script.sh\n']; disp(SBATCH_CMD); JOB2LOAD= JOB2LOAD+1; setup write_sbash_marconi [~,job_info_] = system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); jobid_ = job_info_(21:27); end end end system('rm fort*.90'); disp('done'); -end diff --git a/wk/photomaton.m b/wk/photomaton.m deleted file mode 100644 index 39f6813..0000000 --- a/wk/photomaton.m +++ /dev/null @@ -1,59 +0,0 @@ -if 0 -%% Photomaton : real space - -% Chose the field to plot -% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}'; -% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}' -% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i'; -% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e'; -% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i'; -% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e'; -FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi'; - -% Chose when to plot it -tf = [0 10 50 100]; - -% Slice -ix = 1; iy = 1; iz = 1; -% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')'] -% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')'] -% plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(z(iz))),')'] - -% Averaged -% plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_x'] -% plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_y'] -plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle',FIELDLTX,'\rangle_z'] - - -% -TNAME = []; -fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350]) -plt_2 = @(x) x./max(max(x)); - for i_ = 1:numel(tf) - [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; - subplot(1,numel(tf),i_) - DATA = plt_2(squeeze(plt(FIELD,it))); - pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) - colormap(bluewhitered); caxis([-1,1]); - xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]); - title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); - end - legend(latexize(FIELDLTX)); -save_figure -end - -if 0 -%% Photomaton : quiver ExB velocity -figure -skip = 2; -plt = @(x) x./max(max(x)); -FNAME = 'ZF'; FIELDLTX = '$\bm{u}^{ZF}$'; X_XY = RR; Y_XY = ZZ; -tf = 1200; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)]; -UY = plt(drphi(1:skip:end,1:skip:end,it1)); UX = plt(-dzphi(1:skip:end,1:skip:end,it1)); -pclr = pcolor(X_XY,Y_XY,plt(ni00(:,:,it1))); set(pclr, 'edgecolor','none'); -hold on -quiver((X_XY(1:skip:end,1:skip:end)),(Y_XY(1:skip:end,1:skip:end)),UX,UY,'r'); xlim(L/2*[-1 1]); ylim(L/2*[-1 1]); -pbaspect([1 1 1]) -xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); -title(sprintf('$t c_s/R=%.0f$',Ts2D(it1))); -end \ No newline at end of file diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m index 8700441..a86a524 100644 --- a/wk/shearless_linear_fluxtube.m +++ b/wk/shearless_linear_fluxtube.m @@ -1,110 +1,157 @@ RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.0; % Collision frequency TAU = 1.0; % e/i temperature ratio K_N = 2.22; % Density gradient drive K_T = 6.0; % Temperature ''' SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 1; % Kinetic (1) or adiabatic (0) electron model %% GRID PARAMETERS -NX = 1; % real space x-gridpoints -NY = 2; % '' y-gridpoints -LX = 0; % Size of the squared frequency domain +NX = 2; % real space x-gridpoints +NY = 16; % '' y-gridpoints +LX = 2*pi/0.25; % Size of the squared frequency domain LY = 2*pi/0.25; % Size of the squared frequency domain NZ = 24; % number of perpendicular planes (parallel grid) Q0 = 2.7; % safety factor SHEAR = 0.0; % magnetic shear EPS = 0.18; % inverse aspect ratio +SG = 1; % Staggered z grids option %% TIME PARMETERS -TMAX = 10; % Maximal time unit -DT = 1e-3; % Time step +TMAX = 40; % Maximal time unit +DT = 8e-3; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 10; % 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 JOB2LOAD= -1; %% OPTIONS SIMID = 'shearless_fluxtube'; % Name of the simulation % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK) CO = 1; INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax)) NL_CLOS =-1; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) %% OUTPUTS W_DOUBLE = 0; 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.5; % Hyper diffusivity cutoff ratio -kmax = Nx*pi/Lx;% Highest fourier mode +kmax = NX*pi/LX;% Highest fourier mode NU_HYP = 0.0; % Hyperdiffusivity coefficient MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f K_E = 0.00; % Electrostat ''' GRADB = 1.0; CURVB = 1.0; INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; INIT_PHI= 0; NOISE0 = 0.0e-4; % Init noise amplitude BCKGD0 = 1.0e-4; % Init background LAMBDAD = 0.0; KXEQ0 = 0; % put kx = 0 NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1) %% PARAMETER SCANS -if 1 %% Parameter scan over PJ % PA = [2 4]; % JA = [1 2]; PA = [4]; JA = [2]; DTA= DT*ones(size(JA));%./sqrt(JA); % DTA= DT; mup_ = MU_P; muj_ = MU_J; Nparam = numel(PA); param_name = 'PJ'; -gamma_Ni00 = zeros(Nparam,floor(Nx/2)+1); -gamma_Nipj = zeros(Nparam,floor(Nx/2)+1); -gamma_phi = zeros(Nparam,floor(Nx/2)+1); -% Ni00_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D); -% PHI_ST = zeros(Nparam,floor(Nx/2)+1,TMAX/SPS3D); +Nkx = NX; +Nky = NY/2+1; +gamma_Ni00 = zeros(Nparam,Nkx,Nky); +gamma_Nipj = zeros(Nparam,Nkx,Nky); +gamma_phi = zeros(Nparam,Nkx,Nky); for i = 1:Nparam % Change scan parameter PMAXE = PA(i); PMAXI = PA(i); JMAXE = JA(i); JMAXI = JA(i); DT = DTA(i); setup % system(['rm fort*.90']); % Run linear simulation - if RUN + if 0 system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk']) end % Load and process results %% filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results + + for ikx = 1:Nkx + for iky = 1:Nky + tend = max(Ts3D(abs(Ni00(ikx,iky,1,:))~=0)); + tstart = 0.2*tend; + [~,itstart] = min(abs(Ts3D-tstart)); + [~,itend] = min(abs(Ts3D-tend)); + trange = itstart:itend; + % exp fit on moment 00 + X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,iky,1,trange))); + gamma_Ni00(i,ikx,iky) = LinearFit_s(X_,Y_); + % exp fit on phi + X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,iky,1,trange))); + gamma_phi (i,ikx,iky) = LinearFit_s(X_,Y_); + end + end + gamma_Ni00(i,:,:) = real(gamma_Ni00(i,:,:));% .* (gamma_Ni00(i,:)>=0.0)); + gamma_Nipj(i,:,:) = real(gamma_Nipj(i,:,:));% .* (gamma_Nipj(i,:)>=0.0)); + if 0 + %% Fit verification + figure; + for i = 1:1:Nky + X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(1,i,1,:))); + semilogy(X_,Y_,'DisplayName',['k=',num2str(ky(i))]); hold on; + end + end end +if 1 +%% Plot +SCALE = 1;%sqrt(2); +fig = figure; FIGNAME = 'linear_study'; +plt = @(x) squeeze(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*ky(1:Nky)),plt(gamma_phi(i,2,1:end)),... + 'Color',clr,... + 'LineStyle',linestyle{1},'Marker','^',... +...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']); + 'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']); + hold on; + end + grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(ky)]); + title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$']) + legend('show'); %xlim([0.01,10]) +saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']); +saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']); end if 0 %% Trajectories of some modes figure; -for i = 1:10:Nx/2+1 - semilogy(Ts3D,squeeze(abs(Ne00(i,2,1,:))),'DisplayName',['k=',num2str(kx(i))]); hold on; +for i = 1:1:Nky + semilogy(Ts3D,squeeze(abs(Ni00(1,i,1,:))),'DisplayName',['k=',num2str(ky(i))]); hold on; end end \ No newline at end of file