diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m index dea4c3f..f9ddf90 100644 --- a/matlab/plot/spectrum_1D.m +++ b/matlab/plot/spectrum_1D.m @@ -1,141 +1,141 @@ function [ FIGURE ] = spectrum_1D( data, options ) options.PLAN = 'kxky'; options.COMP = 'avg'; options.INTERP = 0; options.POLARPLOT = 0; options.AXISEQUAL = 1; toplot = process_field(data,options); t = data.Ts3D; frames = toplot.FRAMES; colors = jet(numel(frames)); switch options.NAME case '\Gamma_x' FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS]; yname = '$\sum_{k_y}|\Gamma_k|$'; fieldname = 'transport'; case '\phi' FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS]; yname = '$\sum_{k_y}|\phi|$'; fieldname = 'ES pot.'; case 'n_i' FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS]; yname = '$\sum_{k_y}|n_i|$'; fieldname = 'ion dens.'; case 'n_e' FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS]; yname = '$\sum_{k_y}|n_e|$'; fieldname = 'elec. dens.'; end PLOT2D = 0; switch options.COMPXY case 'avg' - compx = @(x) mean(x,1); - compy = @(x) mean(x,2); + compx = @(x) mean(x,2); + compy = @(x) mean(x,1); case 'sum' - compx = @(x) sum(x,1); - compy = @(x) sum(x,2); + compx = @(x) sum(x,2); + compy = @(x) sum(x,1); case 'max' - compx = @(x) max(x,1); - compy = @(x) max(x,2); + compx = @(x) max(x,2); + compy = @(x) max(x,1); otherwise compx = @(x) x(:,:); compy = @(x) x(:,:); PLOT2D= 1; end if ~PLOT2D - set(gcf, 'Position', [20 50 5000 2000]) + set(gcf, 'Position', [20 50 1200 500]) subplot(1,2,1) - k = data.kx; - xname = '$k_x$'; + k = data.ky; + xname = '$k_y$'; nmax = ceil(data.Nkx*2/3); shiftx = @(x) x;%(1:nmax); shifty = @(x) x;%(1:nmax); switch options.COMPT case 'avg' it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); - Gk = compy(abs(mean(toplot.FIELD(:,:,it0:it1),3))); + Gk = compx(abs(mean(toplot.FIELD(:,:,:),3))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = shiftx(k); if options.OK Y = shifty(Gk)./X; else Y = shifty(Gk); end plot(X,Y,'DisplayName','t-averaged') otherwise for it = 1:numel(toplot.FRAMES) Gk = compy(abs(toplot.FIELD(:,:,toplot.FRAMES(it)))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = shiftx(k); if options.OK Y = shifty(Gk)./X; else Y = shifty(Gk); end plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],... 'Color',colors(it,:)); hold on; end end grid on title(['HeLaZ $k_x$ ',fieldname,' spectrum']); legend('show','Location','eastoutside') xlabel(xname); ylabel(yname) subplot(1,2,2) - k = data.ky; - xname = '$k_y$'; + k = data.kx; + xname = '$k_x$'; nmax = floor(data.Nky/2*2/3); switch options.COMPT case 'avg' - it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); - Gk = compx(abs(mean(toplot.FIELD(:,:,it0:it1),3)))'; +% it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); + Gk = compy(abs(mean(toplot.FIELD(:,:,:),3)))'; Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = k(k>0); X = X(1:end-1); Yp= Gk(k>0); Ym= Gk(k<0); Y = Yp(1:end-1) + Ym(end:-1:1); Y = Y(end:-1:1); plot(X,Y,'DisplayName','t-averaged') otherwise for it = 1:numel(toplot.FRAMES) - Gk = compx(abs(toplot.FIELD(:,:,toplot.FRAMES(it)))); + Gk = compx(abs(toplot.FIELD(:,:,it))); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end X = k(k>0); X = X(1:end-1); Yp= Gk(k>0); Ym= Gk(k<0); Y = Yp(1:end-1) + Ym(end:-1:1); Y = Y(end:-1:1); plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],... 'Color',colors(it,:)); hold on; end end grid on % title('HeLaZ $k_y$ transport spectrum'); legend('show','Location','eastoutside'); xlabel(xname); ylabel(yname) else - it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); - Gk = mean(abs(toplot.FIELD(:,:,it0:it1)),3); +% it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end); + Gk = mean(abs(toplot.FIELD(:,:,:)),3); Gk = squeeze(Gk); if options.NORM Gk = Gk./max(max(abs(Gk))); end pclr = pcolor(toplot.X,toplot.Y,Gk); set(pclr, 'edgecolor','none'); end diff --git a/src/diagnose.F90 b/src/diagnose.F90 index 894263d..c1b81f2 100644 --- a/src/diagnose.F90 +++ b/src/diagnose.F90 @@ -1,116 +1,116 @@ SUBROUTINE diagnose(kstep) ! Diagnostics, writing simulation state to disk USE basic USE diagnostics_par IMPLICIT NONE INTEGER, INTENT(in) :: kstep CALL cpu_time(t0_diag) ! Measuring time !! Basic diagnose loop for reading input file, displaying advancement and ending IF ((kstep .EQ. 0)) THEN INQUIRE(unit=lu_in, name=input_fname) CLOSE(lu_in) ENDIF IF (kstep .GE. 0) THEN ! Terminal info IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax ENDIF ELSEIF (kstep .EQ. -1) THEN CALL cpu_time(finish) ! Display computational time cost IF (my_id .EQ. 0) CALL display_h_min_s(finish-start) END IF !! Specific diagnostic calls - ! CALL diagnose_full(kstep) - IF(nsave_5d .GT. 0) CALL diagnose_moments(kstep) - IF(nsave_3d .GT. 0) CALL diagnose_momspectrum(kstep) - IF(nsave_3d .GT. 0) CALL diagnose_fields(kstep) - IF(nsave_0d .GT. 0) CALL diagnose_profiler(kstep) - IF(nsave_0d .GT. 0) CALL diagnose_gridgeom(kstep) - IF(nsave_0d .GT. 0) CALL diagnose_timetraces(kstep) + CALL diagnose_full(kstep) + ! IF(nsave_5d .GT. 0) CALL diagnose_moments(kstep) + ! IF(nsave_3d .GT. 0) CALL diagnose_momspectrum(kstep) + ! IF(nsave_3d .GT. 0) CALL diagnose_fields(kstep) + ! IF(nsave_0d .GT. 0) CALL diagnose_profiler(kstep) + ! IF(nsave_0d .GT. 0) CALL diagnose_gridgeom(kstep) + ! IF(nsave_0d .GT. 0) CALL diagnose_timetraces(kstep) CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) END SUBROUTINE diagnose SUBROUTINE init_outfile(comm,file0,file,fid) USE diagnostics_par, ONLY : write_doubleprecision, diag_par_outputinputs, input_fname USE basic, ONLY : my_id, jobnum, lu_in, basic_outputinputs USE grid, ONLY : grid_outputinputs USE geometry, ONLY : geometry_outputinputs USE model, ONLY : model_outputinputs USE collision, ONLY : coll_outputinputs USE initial_par, ONLY : initial_outputinputs USE time_integration,ONLY : time_integration_outputinputs USE futils, ONLY : creatf, creatg, creatd, attach, putfile IMPLICIT NONE !input INTEGER, INTENT(IN) :: comm CHARACTER(len=256), INTENT(IN) :: file0 CHARACTER(len=256), INTENT(OUT) :: file INTEGER, INTENT(OUT) :: fid CHARACTER(len=256) :: str,fname INCLUDE 'srcinfo.h' ! Writing output filename WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0) ,'_',jobnum,'.h5' ! 1.1 Initial run ! Main output file creation IF (write_doubleprecision) THEN CALL creatf(file, fid, real_prec='d', mpicomm=comm) ELSE CALL creatf(file, fid, mpicomm=comm) END IF IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(file), ' created' ! basic data group CALL creatg(fid, "/data", "data") ! File group CALL creatg(fid, "/files", "files") CALL attach(fid, "/files", "jobnum", jobnum) ! Add input namelist variables as attributes of /data/input, defined in srcinfo.h ! IF (my_id .EQ. 0) WRITE(*,*) 'VERSION=', VERSION ! IF (my_id .EQ. 0) WRITE(*,*) 'BRANCH=', BRANCH ! IF (my_id .EQ. 0) WRITE(*,*) 'AUTHOR=', AUTHOR ! IF (my_id .EQ. 0) WRITE(*,*) 'HOST=', HOST ! Add the code info and parameters to the file WRITE(str,'(a,i2.2)') "/data/input" CALL creatd(fid, 0,(/0/),TRIM(str),'Input parameters') CALL attach(fid, TRIM(str), "version", VERSION) !defined in srcinfo.h CALL attach(fid, TRIM(str), "branch", BRANCH) !defined in srcinfo.h CALL attach(fid, TRIM(str), "author", AUTHOR) !defined in srcinfo.h CALL attach(fid, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h CALL attach(fid, TRIM(str), "host", HOST) !defined in srcinfo.h CALL basic_outputinputs(fid,str) CALL grid_outputinputs(fid, str) CALL geometry_outputinputs(fid, str) CALL diag_par_outputinputs(fid, str) CALL model_outputinputs(fid, str) CALL coll_outputinputs(fid, str) CALL initial_outputinputs(fid, str) CALL time_integration_outputinputs(fid, str) ! Save STDIN (input file) of this run IF(jobnum .LE. 99) THEN WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum ELSE WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum END IF CALL putfile(fid, TRIM(str), TRIM(input_fname),ionode=0) END SUBROUTINE init_outfile !! Auxiliary routines hidden in headers -! INCLUDE 'diag_headers/diagnose_full.h' -INCLUDE 'diag_headers/diagnose_moments.h' -INCLUDE 'diag_headers/diagnose_momspectrum.h' -INCLUDE 'diag_headers/diagnose_fields.h' -INCLUDE 'diag_headers/diagnose_profiler.h' -INCLUDE 'diag_headers/diagnose_gridgeom.h' -INCLUDE 'diag_headers/diagnose_timetraces.h' +INCLUDE 'diag_headers/diagnose_full.h' +! INCLUDE 'diag_headers/diagnose_moments.h' +! INCLUDE 'diag_headers/diagnose_momspectrum.h' +! INCLUDE 'diag_headers/diagnose_fields.h' +! INCLUDE 'diag_headers/diagnose_profiler.h' +! INCLUDE 'diag_headers/diagnose_gridgeom.h' +! INCLUDE 'diag_headers/diagnose_timetraces.h' diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index 07f4a3d..0b8fef0 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,204 +1,204 @@ addpath(genpath([helazdir,'matlab'])) % ... add addpath(genpath([helazdir,'matlab/plot'])) % ... add addpath(genpath([helazdir,'matlab/compute'])) % ... add addpath(genpath([helazdir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [helazdir,'results/',outfile,'/']; MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; system(['mkdir -p ',MISCDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 10; data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) options.TAVG_0 = 0.8*data.Ts3D(end); options.TAVG_1 = data.Ts3D(end); % Averaging times duration options.NMVA = 1; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 1; fig = plot_radial_transport_and_spacetime(data,options); save_figure(data,fig) end if 0 %% statistical transport averaging options.T = [16000 17000]; fig = statistical_transport_averaging(data,options); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 0; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 9; % options.TIME = dat.Ts5D; -options.TIME = 300:1:4500; +options.TIME = 900:1:990; data.EPS = 0.1; data.a = data.EPS * 2000; create_film(data,options,'.gif') end if 0 %% 2D snapshots % Options -options.INTERP = 0; +options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; -options.NAME = '\phi'; +% options.NAME = '\phi'; % options.NAME = 'n_i'; -% options.NAME = 'N_i^{00}'; +options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; -options.PLAN = 'kxky'; +options.PLAN = 'xy'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; -options.TIME = [300 350 400]; +options.TIME = [900 923 927 990]; data.a = data.EPS * 2e3; fig = photomaton(data,options); save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 1; options.NAME = 'n_i'; options.PLANES = [1:2:12]; options.TIME = [60]; options.PLT_MTOPO = 1; data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig) end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) % options.SPAR = linspace(-3,3,64)+(6/127/2); % options.XPERP = linspace( 0,6,64); options.SPAR = gene_data.vp'; options.XPERP = gene_data.mu'; options.Z = 'avg'; options.T = 900; options.PLT_FCT = 'contour'; options.ONED = 0; options.non_adiab = 1; options.SPECIE = 'i'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); save_figure(data,fig) end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; -options.ST = 0; +options.ST = 1; options.PLOT_TYPE = 'space-time'; -options.NORMALIZED = 0; +options.NORMALIZED = 1; options.JOBNUM = 0; -options.TIME = [300 500]; +options.TIME = [900 950]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% Time averaged spectrum -options.TIME = 1000:5000; -options.NORM = 0; +options.TIME = 1000:1200; +options.NORM =1; options.NAME = '\phi'; -options.NAME = 'n_i'; -options.NAME ='\Gamma_x'; +% options.NAME = 'n_i'; +% options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; -options.COMPXY = 0; +options.COMPXY = 'avg'; options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig) end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 'avg'; options.COMPY = 'avg'; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 1; options.K2PLOT = 1; options.TIME = 5:1:15; options.NMA = 1; options.NMODES = 5; options.iz = 1; fig = mode_growth_meter(data,options); save_figure(data,fig) end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 200; TAVG_1 = 3000; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig) end if 0 %% Metric infos fig = plot_metric(data); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 15]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end diff --git a/wk/analysis_header.m b/wk/analysis_header.m index 012bd79..a97d0de 100644 --- a/wk/analysis_header.m +++ b/wk/analysis_header.m @@ -1,16 +1,16 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" helazdir = '/home/ahoffman/HeLaZ/'; % Directory of the simulation (from results) % if 1% Local results outfile =''; outfile =''; outfile =''; % outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0'; -outfile ='shearless_cyclone/128x128x16xdmax_8_L_120_CBC_1.0'; +outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0'; % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK'; % outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1'; % outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK'; % outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/'; % outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/'; % outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1'; run analysis_3D