diff --git a/matlab/setup.m b/matlab/setup.m index cf9e1f1..71b95cd 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,190 +1,190 @@ %% _______________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nx = NX; % x grid resolution GRID.Lx = LX; % x length GRID.Nexc = NEXC; % to extend Lx when s>0 GRID.Ny = NY; % y '' GRID.Ly = LY; % y '' GRID.Nz = NZ; % z resolution GRID.Npol = NPOL; % z resolution if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end; % Geometry GEOM.geom = ['''',GEOMETRY,'''']; GEOM.q0 = Q0; % q factor GEOM.shear = SHEAR; % shear GEOM.eps = EPS; % inverse aspect ratio GEOM.kappa = KAPPA; % elongation GEOM.delta = DELTA; % triangularity GEOM.zeta = ZETA; % squareness % Model parameters MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; MODEL.LINEARITY = ['''',LINEARITY,'''']; MODEL.KIN_E = KIN_E; if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end; MODEL.beta = BETA; MODEL.mu_x = MU_X; MODEL.mu_y = MU_Y; MODEL.N_HD = N_HD; MODEL.mu_z = MU_Z; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % mass ratio sqrt(m_a/m_i) MODEL.sigma_e = SIGMA_E; MODEL.sigma_i = 1.0; % charge q_a/e MODEL.q_e =-1.0; MODEL.q_i = 1.0; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.K_Ni = K_Ni; MODEL.K_Ne = K_Ne; MODEL.K_Ti = K_Ti; MODEL.K_Te = K_Te; MODEL.GradB = GRADB; % Magnetic gradient MODEL.CurvB = CURVB; % Magnetic curvature MODEL.lambdaD = LAMBDAD; % Collision parameters COLL.collision_model = ['''',CO,'''']; if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end; if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end; COLL.mat_file = '''null'''; switch CO case 'SG' COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5'''; case 'LR' COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''; case 'LD' COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''; % COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5'''; % COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5'''; % COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5'''; end COLL.coll_kcut = COLL_KCUT; % Time integration and intialization parameters -TIME_INTEGRATION.numerical_scheme = '''RK4'''; +TIME_INTEGRATION.numerical_scheme = ['''',NUMERICAL_SCHEME,'''']; INITIAL.INIT_OPT = ['''',INIT_OPT,'''']; INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,'''']; INITIAL.init_background = BCKGD0; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; % Naming and creating input file CONAME = CO; if GKCO CONAME = [CONAME,'GK']; else CONAME = [CONAME,'DK']; end if ~ABCO CONAME = [CONAME,'aa']; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end % Hermite-Laguerre degrees naming if (PMAXE == PMAXI) && (JMAXE == JMAXI) HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)]; else HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),... '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)]; end % temp. dens. drives drives_ = []; if abs(K_Ni) > 0; drives_ = [drives_,'_kN_',num2str(K_Ni)]; end; if abs(K_Ti) > 0; drives_ = [drives_,'_kT_',num2str(K_Ti)]; end; % collision coll_ = ['_nu_%1.1e_',CONAME]; coll_ = sprintf(coll_,NU); % nonlinear lin_ = []; if ~LINEARITY; lin_ = '_lin'; end adiabe_ = []; if ~KIN_E; adiabe_ = '_adiabe'; end % resolution and boxsize res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)]; if (LX ~= LY) geo_ = ['_Lx_',num2str(LX),'_Ly_',num2str(LY)]; else geo_ = ['_L_',num2str(LX)]; end if (NZ > 1) %3D case res_ = [res_,'x',num2str(NZ)]; if abs(Q0) > 0 geo_ = [geo_,'_q0_',num2str(Q0)]; end if abs(EPS) > 0 geo_ = [geo_,'_e_',num2str(EPS)]; end if abs(SHEAR) > 0 geo_ = [geo_,'_s_',num2str(SHEAR)]; end end switch GEOMETRY case 'circular' geo_ = [geo_,'_circ_']; end % put everything together in the param character chain u_ = '_'; % underscore variable PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; if W_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.false.';end; if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end; if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.false.';end; OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) end if ~exist(BASIC.RESDIR, 'dir') mkdir(BASIC.RESDIR) end if ~exist(BASIC.MISCDIR, 'dir') mkdir(BASIC.MISCDIR) end %% Compile and WRITE input file INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; % system(MAKE); %% disp(['Set up ',SIMID]); disp([res_,geo_,HLdeg_]); if JOB2LOAD>=0 disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m index e637d9d..b7bffaf 100644 --- a/wk/analysis_gyacomo.m +++ b/wk/analysis_gyacomo.m @@ -1,233 +1,233 @@ %% UNCOMMENT FOR TUTORIAL gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory % resdir = '.'; %Name of the directory where the results are located % JOBNUMMIN = 00; JOBNUMMAX = 10; %% addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add addpath(genpath([gyacomodir,'matlab/compute'])) % ... add addpath(genpath([gyacomodir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [gyacomodir,resdir,'/']; MISCDIR = ['/misc/',resdir,'/']; %For long term storage system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) % data.scale = 1;%/(data.Nx*data.Ny)^2; i_ = 1; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.TAVG_0 = data.TJOB_SE(i_)+600;%0.4*data.Ts3D(end); options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration options.NCUT = 4; % Number of cuts for averaging and error estimation options.NMVA = 100; % 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 = 0; options.RESOLUTION = 256; fig = plot_radial_transport_and_spacetime(data,options); % save_figure(data,fig,'.png') end if 0 %% statistical transport averaging gavg =[]; gstd = []; for i_ = 3:2:numel(data.TJOB_SE) % i_ = 3; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)]; options.NPLOTS = 0; [fig, gavg_, gstd_] = statistical_transport_averaging(data,options); gavg = [gavg gavg_]; gstd = [gstd gstd_]; end disp(gavg); disp(gstd); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; % options.NAME = '\phi'; % options.NAME = '\omega_z'; -options.NAME = 'N_e^{00}'; +options.NAME = 'N_i^{00}'; % options.NAME = 'v_x'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; % options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts3D; options.TIME = [000:1:8000]; data.EPS = 0.1; data.a = data.EPS * 2000; options.RESOLUTION = 256; create_film(data,options,'.gif') end if 0 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = '\psi'; % options.NAME = 'n_i'; % options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; -options.PLAN = 'xy'; +options.PLAN = 'kxky'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [1000 3000 5000 8000 10000]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 0; options.NAME = '\phi'; options.PLANES = [1]; options.TIME = [30]; options.PLT_MTOPO = 1; options.PLT_FTUBE = 0; data.EPS = 0.4; data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) % options.SPAR = linspace(-3,3,32)+(6/127/2); % options.XPERP = linspace( 0,6,32); options.SPAR = gene_data.vp'; options.XPERP = gene_data.mu'; options.iz = 'avg'; options.T = [1500 5000]; options.PLT_FCT = 'contour'; options.ONED = 1; options.non_adiab = 0; options.SPECIE = 'e'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [1000]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); % save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum options.TIME = [2000 3000]; options.NORM =1; % options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 1; options.COMPY = 2; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution -options.NORMALIZED = 0; -options.K2PLOT = [0.1 0.2 0.3 0.4]; +options.NORMALIZED = 1; +options.K2PLOT = [0.1 0.5 1.0 2.0 6.0]; options.TIME = [00:1200]; options.NMA = 1; options.NMODES = 5; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos options.SHOW_FLUXSURF = 1; options.SHOW_METRICS = 0; fig = plot_metric(data,options); 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 20]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); % save_figure(data,fig,'.png') end diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index 5e9a28c..aa7ca49 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,71 +1,71 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" gyacomodir = '/home/ahoffman/gyacomo/'; % Directory of the simulation (from results) % if 1% Local results % resdir ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits % resdir ='shearless_cyclone/128x64x16x5x3_Dim_90'; % resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; % resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; % resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; %% shearless CBC % resdir ='shearless_cyclone/64x32x16x5x3_CBC_080'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_120'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_080'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_100'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_120'; % resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% CBC % resdir = 'CBC/64x32x16x5x3'; % resdir = 'CBC/64x128x16x5x3'; % resdir = 'CBC/128x64x16x5x3'; % resdir = 'CBC/96x96x16x3x2_Nexc_6'; % resdir = 'CBC/128x96x16x3x2'; % resdir = 'CBC/192x96x16x3x2'; % resdir = 'CBC/192x96x24x13x7'; % resdir = 'CBC/kT_11_128x64x16x5x3'; % resdir = 'CBC/kT_9_256x128x16x3x2'; % resdir = 'CBC/kT_4.5_128x64x16x13x3'; % resdir = 'CBC/kT_4.5_192x96x24x13x7'; % resdir = 'CBC/kT_4.5_128x64x16x13x7'; % resdir = 'CBC/kT_4.5_128x96x24x15x5'; % resdir = 'CBC/kT_5.3_192x96x24x13x7'; % resdir = 'CBC/kT_13_large_box_128x64x16x5x3'; % resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02'; % resdir = 'CBC/kT_scan_128x64x16x5x3'; % resdir = 'CBC/kT_scan_192x96x16x3x2'; % resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6'; % resdir = 'dbg/nexc_dbg'; % resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4'; % resdir = 'CBC_Ke_EM/192x96x24x5x3'; % resdir = 'CBC_Ke_EM/96x48x16x5x3'; % resdir = 'CBC_Ke_EM/minimal_res'; %% KBM % resdir = 'NL_KBM/192x64x24x5x3'; %% Linear CBC % resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; % resdir = 'testcases/miller_example'; % resdir = 'Miller/128x256x3x2_CBC_dt_5e-3'; %% CBC Miller -% resdir = 'GCM_CBC/daint/Miller_GX_comparison'; +resdir = 'GCM_CBC/daint/Miller_GX_comparison'; %% RK3 -resdir = 'dbg/SSPx_RK3_test'; -resdir = 'dbg/SSPx_RK3_test/RK4'; -resdir = ['results/',resdir]; +% resdir = 'dbg/SSPx_RK3_test'; +% resdir = 'dbg/SSPx_RK3_test/RK4'; +% resdir = ['results/',resdir]; JOBNUMMIN = 00; JOBNUMMAX = 10; run analysis_gyacomo diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m index 9b6c383..cdc2a5f 100644 --- a/wk/lin_ETPY.m +++ b/wk/lin_ETPY.m @@ -1,187 +1,187 @@ %% QUICK RUN SCRIPT for linear entropy mode (ETPY) in a Zpinch % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % -% SIMID = 'lin_ETPY'; % Name of the simulation -SIMID = 'dbg'; % Name of the simulation +SIMID = 'lin_ETPY'; % Name of the simulation +% SIMID = 'dbg'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options PROGDIR = '/home/ahoffman/gyacomo/'; -% EXECNAME = 'gyacomo_dbg'; -EXECNAME = 'gyacomo'; +EXECNAME = 'gyacomo_dbg'; +% EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.05; % Collision frequency +NU = 0.01; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 2.0; % ele Density ''' -K_Te = 0.4; % ele Temperature ''' -K_Ni = 2.0; % ion Density gradient drive -K_Ti = 0.4; % ion Temperature ''' +K_Ne = 2.5; % ele Density ''' +K_Te = K_Ne/4; % ele Temperature ''' +K_Ni = K_Ne; % ion Density gradient drive +K_Ti = K_Ni/4; % ion Temperature ''' SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS -P = 4; +P = 8; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 2; % real space x-gridpoints -NY = 50; % '' y-gridpoints +NY = 100; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 2*pi/0.05; % Size of the squared frequency domain +LY = 200;%2*pi/0.05; % Size of the squared frequency domain NZ = 1; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -GEOMETRY= 'Z-pinch'; +GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear -KAPPA = 1.0; % elongation +KAPPA = 0.0; % elongation DELTA = 0.0; % triangularity ZETA = 0.0; % squareness NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 100; % Maximal time unit +TMAX = 500; % 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 +SPS2D = -1; % Sampling per time unit for 2D arrays SPS3D = 1; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; -GKCO = 0; % gyrokinetic operator +GKCO = 1; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) -INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom +INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom +NUMERICAL_SCHEME = 'RK2'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5 %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; -MU_Z = 0.0; % +MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 1.0e-5; % Init noise amplitude BCKGD0 = 0.0; % Init background GRADB = 1.0; CURVB = 1.0; COLL_KCUT = 1000; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [PROGDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) options.TRANGE = [0.5 1]*data.Ts3D(end); options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3 lg = compute_fluxtube_growth_rate(data,options); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 0 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.9]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 10; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 10; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m index ebca810..339b183 100644 --- a/wk/lin_ITG.m +++ b/wk/lin_ITG.m @@ -1,191 +1,193 @@ %% QUICK RUN SCRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % % SIMID = 'lin_ITG'; % Name of the simulation SIMID = 'dbg'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options gyacomodir = '/home/ahoffman/gyacomo/'; % EXECNAME = 'gyacomo_1.0'; % EXECNAME = 'gyacomo_dbg'; EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.05; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 0*2.22; % ele Density ''' +K_Ne = 2.22; % ele Density ''' K_Te = 6.96; % ele Temperature ''' -K_Ni = 0*2.22; % ion Density gradient drive +K_Ni = 2.22; % ion Density gradient drive K_Ti = 6.96; % ion Temperature ''' SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 20; % real space x-gridpoints -NY = 2; % '' y-gridpoints +NY = 20; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 2*pi/0.3; % Size of the squared frequency domain +LY = 2*pi/0.05; % Size of the squared frequency domain NZ = 32; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -% GEOMETRY= 's-alpha'; -GEOMETRY= 'miller'; +GEOMETRY= 's-alpha'; +% GEOMETRY= 'miller'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear KAPPA = 1.0; % elongation DELTA = 0.0; % triangularity ZETA = 0.0; % squareness NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 25; % Maximal time unit -DT = 3e-3; % Time step +TMAX = 20; % 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 +SPS2D = -1; % Sampling per time unit for 2D arrays SPS3D = 5; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom +NUMERICAL_SCHEME = 'RK3'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5 %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; -NOISE0 = 0.0e-5; % Init noise amplitude -BCKGD0 = 1.0; % Init background +NOISE0 = 1.0e-5; % Init noise amplitude +BCKGD0 = 0.0; % Init background GRADB = 1.0; CURVB = 1.0; +COLL_KCUT = 1000; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN % system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [gyacomodir,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) options.TRANGE = [0.5 1]*data.Ts3D(end); options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3 lg = compute_fluxtube_growth_rate(data,options); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 1 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.3]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 1 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end diff --git a/wk/lin_RHT.m b/wk/lin_RHT.m index dcdcc75..d9a2300 100644 --- a/wk/lin_RHT.m +++ b/wk/lin_RHT.m @@ -1,186 +1,189 @@ %% QUICK RUN CRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % SIMID = 'lin_RHT'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options -HELAZDIR = '/home/ahoffman/HeLaZ/'; -% EXECNAME = 'helaz3'; -EXECNAME = 'helaz3'; +GYACOMODIR = '/home/ahoffman/gyacomo/'; +% EXECNAME = 'gyacomo_dbg'; +EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.0; % Collision frequency +NU = 0.1; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 1.0; % ele Density ''' -K_Te = 4.0; % ele Temperature ''' -K_Ni = 2.0; % ion Density gradient drive -K_Ti = 8.0; % ion Temperature ''' +K_Ne = 0.0; % ele Density ''' +K_Te = 0.0; % ele Temperature ''' +K_Ni = 0.0; % ion Density gradient drive +K_Ti = 0.0; % ion Temperature ''' SIGMA_E = 0.5;%0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) -KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons -BETA = 0.05; % electron plasma beta +KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons +BETA = 0.00; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 2; % real space x-gridpoints NY = 2; % '' y-gridpoints -LX = 2*pi/0.0628; % Size of the squared frequency domain -LY = 2*pi/0.1; % Size of the squared frequency domain +LX = pi/2.5; % Size of the squared frequency domain +LY = 2*pi/0.5; % Size of the squared frequency domain NZ = 16; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY % GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps -GEOMETRY= 's-alpha'; -% GEOMETRY= 'circular'; +% GEOMETRY= 's-alpha'; +GEOMETRY= 'miller'; Q0 = 1.4; % safety factor -SHEAR = 0.0; % magnetic shear -NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) +SHEAR = 0.8; % magnetic shear +KAPPA = 1.0; % elongation +DELTA = 0.0; % triangularity +ZETA = 0.0; % squareness +NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) %0: adaptive EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 10; % Maximal time unit -DT = 5e-3; % Time step +TMAX = 50; % 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 = 10; % Sampling per time unit for 2D arrays +SPS2D = -1; % Sampling per time unit for 2D arrays +SPS3D = 5; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays -SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom +NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5 %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 0.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; +COLL_KCUT = 1000; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN -% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) -% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) +% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',GYACOMODIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; -LOCALDIR = [HELAZDIR,'results/',filename,'/']; +LOCALDIR = [GYACOMODIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 0 %% linear growth rate (adapted for 2D zpinch and fluxtube) trange = [0.5 1]*data.Ts3D(end); nplots = 3; lg = compute_fluxtube_growth_rate(data,trange,nplots); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 0 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.1]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; -options.K2PLOT = 1; -options.TIME = [0:1000]; +options.K2PLOT = 2; +options.TIME = [0:50]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); -save_figure(data,fig,'.png') +% save_figure(data,fig,'.png') end if 1 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3));%./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end \ No newline at end of file