diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m index 0eb275e6621f04545cda7df520fc2f4520732f0d..f549d7b8fed3daaf318e12b01da0b91b120db6ae 100644 --- a/wk/lin_run_script.m +++ b/wk/lin_run_script.m @@ -21,10 +21,22 @@ EXECNAME = 'gyacomo23_sp'; % single precision % EXECNAME = 'gyacomo23_dp'; % double precision %% Setup parameters -% run lin_DTT_AB_rho85_PT +% run lin_DTT_AB_rho85 +% run lin_DTT_AB_rho98 +run lin_JET_rho97 % run lin_Entropy % run lin_ITG - +%% Change parameters +NY = 2; +PMAX = 4; +JMAX = 2; +ky = 10.0; +LY = 2*pi/ky; +DT = 1e-3/ky; +SIGMA_E = 0.023; +TMAX = 1.5/ky; +DTSAVE0D = 0.01; +DTSAVE3D = 0.01; %%------------------------------------------------------------------------- %% RUN setup @@ -32,9 +44,9 @@ setup % Run linear simulation if RUN MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;']; -% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;']; + RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;']; % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;']; - RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 1 4 2 0;']; + % RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 1 4 2 0;']; % RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;']; % RUN = ['./../../../bin/gyacomo23_sp 0;']; MVOUT='cd ../../../wk;'; @@ -53,7 +65,7 @@ data = {}; % Initialize data as an empty cell array data = compile_results_low_mem(data,LOCALDIR,J0,J1); -if 0 % Activate or not +if 1 % Activate or not %% plot mode evolution and growth rates % Load phi [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi'); @@ -70,7 +82,7 @@ options.fftz.flag = 0; % Set fftz.flag option to 0 fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig end -if 1 +if (1 && NZ>4) %% Ballooning plot [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi'); if data.inputs.BETA > 0 diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m index a10d35fa9c82b68f585c143c1b5e0a19e643cd3b..a587fc23207dbe35794a5f3d5bc7092aa80b128d 100644 --- a/wk/lin_scan_script.m +++ b/wk/lin_scan_script.m @@ -16,7 +16,7 @@ addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder %% Setup run or load an executable RUN = 1; % To run or just to load -RERUN = 0; % rerun if the data does not exist +RERUN = 1; % rerun if the data does not exist default_plots_options EXECNAME = 'gyacomo23_sp'; % single precision % EXECNAME = 'gyacomo23_dp'; % double precision @@ -29,13 +29,13 @@ run lin_DTT_AB_rho85_PT %% Modify parameters % NZ = 1; NY = 2; -DT = 2e-3; - +DT = 0.5e-3; +TMAX = 50; +MU_X = 0.1; MU_Y = 0.1; %% Scan parameters SIMID = [SIMID,'_scan']; -P_a = [2 4 8]; -ky_a= 0.1:0.1:0.5; - +P_a = [2 4 6 8]; +ky_a = 0.05:0.05:1.5; %% Scan loop % arrays for the result g_ky = zeros(numel(ky_a),numel(P_a),2); @@ -47,63 +47,63 @@ for PMAX = P_a i = 1; for ky = ky_a LY = 2*pi/ky; - %% RUN - setup - % naming - filename = [SIMID,'/',PARAMS,'/']; - LOCALDIR = [gyacomodir,'results/',filename,'/']; - % check if data exist to run if no data - data_ = {}; - try - data_ = compile_results_low_mem(data_,LOCALDIR,00,00); - Ntime = numel(data_.Ts0D); - catch - data_.outfilenames = []; - end - if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10) - % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk']) - end - data_ = compile_results_low_mem(data_,LOCALDIR,00,00); - [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi'); - if numel(data_.Ts3D)>10 - if numel(data_.Ts3D)>5 - % Load results after trying to run + %% RUN + setup + % naming filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [gyacomodir,'results/',filename,'/']; - + % check if data exist to run if no data + data_ = {}; + try + data_ = compile_results_low_mem(data_,LOCALDIR,00,00); + Ntime = numel(data_.Ts0D); + catch + data_.outfilenames = []; + end + if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10) + % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk']) + end data_ = compile_results_low_mem(data_,LOCALDIR,00,00); [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi'); - - % linear growth rate (adapted for 2D zpinch and fluxtube) - options.TRANGE = [0.5 1]*data_.Ts3D(end); - options.NPLOTS = 0; % 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 - - [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window - [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ... - field = 0; - field_t = 0; - for ik = 2:NY/2+1 - field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z - field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only - to_measure = log(field_t(it1:it2)); - tw = double(data_.Ts3D(it1:it2)); - % gr = polyfit(tw,to_measure,1); - gr = fit(tw,to_measure,'poly1'); - err= confint(gr); - g_ky(i,j,ik) = gr.p1; - g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2; - end - [gmax, ikmax] = max(g_ky(i,j,:)); - - msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg); + if numel(data_.Ts3D)>10 + if numel(data_.Ts3D)>5 + % Load results after trying to run + filename = [SIMID,'/',PARAMS,'/']; + LOCALDIR = [gyacomodir,'results/',filename,'/']; + + data_ = compile_results_low_mem(data_,LOCALDIR,00,00); + [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi'); + + % linear growth rate (adapted for 2D zpinch and fluxtube) + options.TRANGE = [0.5 1]*data_.Ts3D(end); + options.NPLOTS = 0; % 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 + + [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window + [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ... + field = 0; + field_t = 0; + for ik = 2:NY/2+1 + field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z + field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only + to_measure = log(field_t(it1:it2)); + tw = double(data_.Ts3D(it1:it2)); + % gr = polyfit(tw,to_measure,1); + gr = fit(tw,to_measure,'poly1'); + err= confint(gr); + g_ky(i,j,ik) = gr.p1; + g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2; + end + [gmax, ikmax] = max(g_ky(i,j,:)); + + msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg); + end end + i = i + 1; end - i = i + 1; -end -j = j + 1; + j = j + 1; end %% take max growth rate among z coordinate @@ -114,8 +114,11 @@ e_ = g_std(:,:,2); if(numel(ky_a)>1 && numel(P_a)>1) pmin = num2str(min(P_a)); pmax = num2str(max(P_a)); kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a)); - filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,... - '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat']; + filename = [num2str(NX),'x',num2str(NZ),... + '_ky_',kymin,'_',kymax,... + '_P_',pmin,'_',pmax,... + '_kN_',num2str(K_Ni),... + '_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat']; metadata.name = filename; metadata.kymin = ky; metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_Ti),', $\kappa_N=$',num2str(K_Ni)]; diff --git a/wk/parameters/lin_DTT_AB_rho85_PT.m b/wk/parameters/lin_DTT_AB_rho85.m similarity index 92% rename from wk/parameters/lin_DTT_AB_rho85_PT.m rename to wk/parameters/lin_DTT_AB_rho85.m index f950920ec2285b2acbbd676c255369a752a94028..dfdf7fa1bace13a36b1e225ef8f5a0375ba59252 100644 --- a/wk/parameters/lin_DTT_AB_rho85_PT.m +++ b/wk/parameters/lin_DTT_AB_rho85.m @@ -14,8 +14,8 @@ dpdx_gn =0.086; SIMID = 'lin_DTT_AB_rho85_PT'; % Name of the simulation %% Set up physical parameters CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss -nu = nu_ei; %(0.00235 in GENE) -TAU = 0.9360; % e/i temperature ratio +NU = 0.1; %(0.00235 or 0.56 in GENE?) +TAU = 0.9360; % e/i temperature ratio K_Ne = 1.33; % ele Density ''' K_Te = 12.0; % ele Temperature ''' K_Ni = 1.33; % ion Density gradient drive @@ -26,12 +26,11 @@ ADIAB_E = (NA==1); % adiabatic electron model BETA = b_gn; % electron plasma beta MHD_PD = 0; %% Set up grid parameters -P = 4; -J = P/2;%P/2; -PMAX = P; % Hermite basis size -JMAX = J; % Laguerre basis size -NX = 16; % real space x-gridpoints -NY = 16; % real space y-gridpoints + +PMAX = 4; % Hermite basis size +JMAX = PMAX/2; % Laguerre basis size +NX = 4; % real space x-gridpoints +NY = 24; % real space y-gridpoints LX = 2*pi/0.1; % Size of the squared frequency domain in x direction LY = 2*pi/0.2; % Size of the squared frequency domain in y direction NZ = 24; % number of perpendicular planes (parallel grid) @@ -55,7 +54,7 @@ SHIFT_Y = 0.0; % Shift in the periodic BC in z NPOL = 1; % Number of poloidal turns PB_PHASE = 0; %% TIME PARAMETERS -TMAX = 15; % Maximal time unit +TMAX = 25; % Maximal time unit DT = 1e-3; % Time step DTSAVE0D = 0.5; % Sampling time for 0D arrays DTSAVE2D = -1; % Sampling time for 2D arrays diff --git a/wk/parameters/lin_DTT_AB_rho98_PT.m b/wk/parameters/lin_DTT_AB_rho98.m similarity index 91% rename from wk/parameters/lin_DTT_AB_rho98_PT.m rename to wk/parameters/lin_DTT_AB_rho98.m index 55352cd11ae29c5743ed4c35b455bc215998412b..12c49a63686d0bd9a563732f8917fae00b93b94c 100644 --- a/wk/parameters/lin_DTT_AB_rho98_PT.m +++ b/wk/parameters/lin_DTT_AB_rho98.m @@ -7,11 +7,11 @@ mref = 2.0; % in proton mass lnLAMBDA = 13; % Coulomb logarithm nuref = 0.45*2.3031e-5*lnLAMBDA*nref*Lref/Tref/Tref; %(0.00235 in GENE) nu_ei = 0.569013; -nu_gn = 0.00235; +nu_gn = 0.235; b_gn = 0.0039; dpdx_gn =0.086; %% Set simulation parameters -SIMID = 'lin_DTT_AB_rho85_PT'; % Name of the simulation +SIMID = 'lin_DTT_AB_rho98_PT'; % Name of the simulation %% Set up physical parameters CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss nu = nu_ei; %(0.00235 in GENE) @@ -20,7 +20,7 @@ K_Ne = 65; % ele Density ''' K_Te = 350; % ele Temperature ''' K_Ni = K_Ne; % ion Density gradient drive K_Ti = 350; % ion Temperature ''' -SIGMA_E = 0.0233380/sqrt(mref); % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) +SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) NA = 2; % number of kinetic species ADIAB_E = (NA==1); % adiabatic electron model BETA = b_gn; % electron plasma beta @@ -30,11 +30,11 @@ P = 4; J = P/2;%P/2; PMAX = P; % Hermite basis size JMAX = J; % Laguerre basis size -NX = 16; % real space x-gridpoints -NY = 16; % real space y-gridpoints +NX = 8; % real space x-gridpoints +NY = 2; % real space y-gridpoints LX = 2*pi/0.1; % Size of the squared frequency domain in x direction -LY = 2*pi/0.2; % Size of the squared frequency domain in y direction -NZ = 24; % number of perpendicular planes (parallel grid) +LY = 2*pi/0.5; % Size of the squared frequency domain in y direction +NZ = 32; % number of perpendicular planes (parallel grid) SG = 0; % Staggered z grids option NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) diff --git a/wk/parameters/lin_Entropy.m b/wk/parameters/lin_Entropy.m index ea0b655db7e4e66f06ff53dea77bfe562c9ac8ce..31c57f712dffb8b3b1066087d95e926458cafc50 100644 --- a/wk/parameters/lin_Entropy.m +++ b/wk/parameters/lin_Entropy.m @@ -14,17 +14,17 @@ NA = 2; % number of kinetic species ADIAB_E = 0; % adiabatic electron model ADIAB_I = 0; % adiabatic ion model BETA = 0.000; % electron plasma beta -MHD_PD = 0; % MHD pressure drift +MHD_PD = 1; % MHD pressure drift %% Set up grid parameters P = 4; J = P/2;%P/2; PMAX = P; % Hermite basis size JMAX = J; % Laguerre basis size NX = 2; % real space x-gridpoints -NY = 12; % real space y-gridpoints +NY = 40; % real space y-gridpoints LX = 2*pi/0.05; % Size of the squared frequency domain in x direction -LY = 2*pi/0.2; % Size of the squared frequency domain in y direction -NZ = 16; % number of perpendicular planes (parallel grid) +LY = 2*pi/0.1; % Size of the squared frequency domain in y direction +NZ = 1; % number of perpendicular planes (parallel grid) SG = 0; % Staggered z grids option NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) %% GEOMETRY @@ -45,8 +45,8 @@ SHIFT_Y = 0.0; % Shift in the periodic BC in z NPOL = 1; % Number of poloidal turns PB_PHASE= 0; %% TIME PARAMETERS -TMAX = 100; % Maximal time unit -DT = 1e-3; % Time step +TMAX = 50; % Maximal time unit +DT = 1e-2; % Time step DTSAVE0D = 1.0; % Sampling time for 0D arrays DTSAVE2D = -1; % Sampling time for 2D arrays DTSAVE3D = 2.0; % Sampling time for 3D arrays diff --git a/wk/parameters/lin_JET_rho97.m b/wk/parameters/lin_JET_rho97.m new file mode 100644 index 0000000000000000000000000000000000000000..5047e578535c447cb450728dbe55acf7513e2753 --- /dev/null +++ b/wk/parameters/lin_JET_rho97.m @@ -0,0 +1,101 @@ +% Parameters found in Parisi et al. 2020 +% Jet shot 92174 +%% Set simulation parameters +SIMID = 'lin_JET_rho97'; % Name of the simulation +%% Set up physical parameters +CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss +nu = 0.83; %(0.00235 in GENE) +TAU = 0.56; % e/i temperature ratio +K_Ne = 10; % ele Density ''' +K_Te = 42; % ele Temperature ''' +K_Ni = 10; % ion Density gradient drive +K_Ti = 11; % ion Temperature ''' +SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) +NA = 2; % number of kinetic species +ADIAB_E = (NA==1); % adiabatic electron model +BETA = 0.0031; % electron plasma beta +MHD_PD = 0; + +%% GEOMETRY +% GEOMETRY= 's-alpha'; +GEOMETRY= 'miller'; +EPS = 0.9753*0.91/2.91; % inverse aspect ratio +Q0 = 5.10; % safety factor +SHEAR = 3.36; % magnetic shear +KAPPA = 1.55; % elongation +S_KAPPA = 0.95; +DELTA = 0.26; % triangularity +S_DELTA = 0.74; +ZETA = 0; % squareness +S_ZETA = 0; +PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected') +SHIFT_Y = 0.0; % Shift in the periodic BC in z +NPOL = 1; % Number of poloidal turns +PB_PHASE = 0; + +%% Set up grid parameters +P = 4; +J = P/2;%P/2; +PMAX = P; % Hermite basis size +JMAX = J; % Laguerre basis size +NX = 8; % real space x-gridpoints +NY = 2; % real space y-gridpoints +LX = 2*pi/0.1; % Size of the squared frequency domain in x direction +LY = 2*pi/0.5; % Size of the squared frequency domain in y direction +NZ = 32; % number of perpendicular planes (parallel grid) +SG = 0; % Staggered z grids option +NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) + +%% TIME PARAMETERS +TMAX = 15; % Maximal time unit +DT = 1e-3; % Time step +DTSAVE0D = 0.5; % Sampling time for 0D arrays +DTSAVE2D = -1; % Sampling time for 2D arrays +DTSAVE3D = 0.5; % Sampling time for 3D arrays +DTSAVE5D = 100; % Sampling time for 5D arrays +JOB2LOAD = -1; % Start a new simulation serie + +%% OPTIONS +LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) +CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) +GKCO = 1; % Gyrokinetic operator +ABCO = 1; % INTERSPECIES collisions +INIT_ZF = 0; % Initialize zero-field quantities +HRCY_CLOS = 'truncation'; % Closure model for higher order moments +DMAX = -1; +NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments +NMAX = 0; +KERN = 0; % Kernel model (0 : GK) +INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom +NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5) + +%% OUTPUTS +W_DOUBLE = 1; % Output flag for double moments +W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy) +W_HF = 1; % Output flag for high-frequency potential energy +W_PHI = 1; % Output flag for potential +W_NA00 = 1; % Output flag for nalpha00 (density of species alpha) +W_DENS = 1; % Output flag for total density +W_TEMP = 1; % Output flag for temperature +W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha) +W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% UNUSED PARAMETERS +% These parameters are usually not to play with in linear runs +MU = 0.0; % Hyperdiffusivity coefficient +MU_X = MU; % Hyperdiffusivity coefficient in x direction +MU_Y = MU; % Hyperdiffusivity coefficient in y direction +N_HD = 4; % Degree of spatial-hyperdiffusivity +MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction +HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model +MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite +MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre +LAMBDAD = 0.0; % Lambda Debye +NOISE0 = 1.0e-5; % Initial noise amplitude +BCKGD0 = 0.0e-5; % Initial background +k_gB = 1.0; % Magnetic gradient strength +k_cB = 1.0; % Magnetic curvature strength +COLL_KCUT = 1; % Cutoff for collision operator +ADIAB_I = 0; % adiabatic ion model \ No newline at end of file