diff --git a/Makefile b/Makefile
index 12c84caa117960b263ff4464fce0dfc1a9567294..9dbf41acf7da3219f96bf63d79c498c27dbd7d3c 100644
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,11 @@
include local/dirs.inc
include local/make.inc
#Different namings depending on the make input
-EXEC = $(BINDIR)/gyacomo #all
-EFST = $(BINDIR)/gyacomo_fast #fast
+EXEC = $(BINDIR)/gyacomo #all
+EFST = $(BINDIR)/gyacomo_fast #fast
EDBG = $(BINDIR)/gyacomo_debug #debug
EALP = $(BINDIR)/gyacomo_alpha #alpha
-
+EGFT = $(BINDIR)/gyacomo_gfort #gfort version
# F90 = mpiifort
F90 = mpif90
# #F90 = ftn #for piz-daint cluster
@@ -28,13 +28,18 @@ fast: F90FLAGS = -fast
fast: $(EFST)
# Debug version with all flags
debug: dirs src/srcinfo.h
-debug: F90FLAGS = -g -traceback -CB -ftrapuv -warn all -debug all
+debug: F90FLAGS = -g -traceback -ftrapuv -warn all -debug all
+# debug: F90FLAGS = -g -traceback -check all -ftrapuv -warn all -debug all
debug: $(EDBG)
# Alpha version, optimized as all but creates another binary
alpha: dirs src/srcinfo.h
alpha: F90FLAGS = -O3 -xHOST
alpha: $(EALP)
-
+# gfortran version, compile with gfortran
+gfort: dirs src/srcinfo.h
+gfort: F90FLAGS = -g -std=legacy -ffree-line-length-0
+gfort: EXTMOD = -J $(MODDIR)
+gfort: $(EGFT)
install: dirs src/srcinfo.h $(EXEC) mvmod
run: all
@@ -93,6 +98,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(EALP): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
+ $(EGFT): $(FOBJ)
+ $(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
$(OBJDIR)/advance_field_mod.o : src/advance_field_mod.F90 \
$(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o \
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 7d6962c6816c0cf24b46fb728001d72a11b44370..27c14463e467bb13ca11f93c5d559a5eb2920fe3 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -7,6 +7,7 @@ DATA.NU_EVOL = []; % evolution of parameter nu between jobs
DATA.CO_EVOL = []; % evolution of CO
DATA.MUx_EVOL = []; % evolution of parameter mu between jobs
DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
+DATA.MUz_EVOL = []; % evolution of parameter mu between jobs
DATA.K_N_EVOL = []; %
DATA.L_EVOL = []; %
DATA.DT_EVOL = []; %
@@ -53,7 +54,7 @@ while(CONTINUE)
try
openable = ~isempty(h5read(filename,'/data/var3d/time'));
catch
- openable = 0
+ openable = 0;
end
if openable
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -235,6 +236,7 @@ while(CONTINUE)
DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
+ DATA.MUz_EVOL = [DATA.MUz_EVOL DATA.MUz DATA.MUz];
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];
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 076058a77e13513ab248ecdd88dbd45adff55642..1fcc187f8d6070ff4eb27493bd01190d8c9317ff 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -4,14 +4,18 @@ NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
d = OPTIONS.fftz.flag; % To add spectral evolution of z (useful for 3d zpinch)
-
-switch OPTIONS.iz
- case 'avg'
- field = squeeze(mean(DATA.PHI,3));
- zstrcomp = 'z-avg';
- otherwise
- field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
- zstrcomp = ['z=',DATA.z(OPTIONS.iz)];
+if numel(size(DATA.PHI)) == 3
+ field = squeeze(DATA.PHI);
+ zstrcomp = 'z=0';
+else
+ switch OPTIONS.iz
+ case 'avg'
+ field = squeeze(mean(DATA.PHI,3));
+ zstrcomp = 'z-avg';
+ otherwise
+ field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
+ zstrcomp = ['z=',DATA.z(OPTIONS.iz)];
+ end
end
FRAMES = zeros(size(OPTIONS.TIME));
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 2d87dfb19bd0e9e3d9d324bed3197134af50ca46..451428f444e97e9a21f94fbe40ddca677fa6a080 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -20,9 +20,17 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
if cmpx
- data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ if(numel(sz) == 3)
+ data(:,:,it) = tmp.real + 1i * tmp.imaginary;
+ else
+ data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
else
- data(:,:,:,it) = tmp;
+ if(numel(sz) == 3)
+ data(:,:,it) = tmp;
+ else
+ data(:,:,:,it) = tmp;
+ end
end
end
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 6d657839c7a70c7c3b031ba752eb987945113119..1f1b32a8e3444066b28b6ceda3f7cf4b3bfabbdc 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -48,6 +48,7 @@ DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
DATA.MUx = h5readatt(filename,'/data/input','mu_x');
DATA.MUy = h5readatt(filename,'/data/input','mu_y');
+DATA.MUz = h5readatt(filename,'/data/input','mu_z');
try
DATA.BETA = h5readatt(filename,'/data/input','beta');
catch
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 3edaeae6bc6542360906b916261eb9822f41d17b..3ca0baa1e15a6a853f16ba44dd7e1ee895777908 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -78,10 +78,15 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
%---------- CBC
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
-fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
+% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_kine.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_kine.txt';
+%----------Convergence nv CBC
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_8_Lv_3_Lw_6.txt';
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_16_Lv_3_Lw_6.txt';
+fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_24_Lv_3_Lw_6.txt';
+
data_ = load([path,fname]);
figure
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index 2e7ccdcfe3fe061092db9f16d464f9b53e1deba5..a844fafe64824c2347266edf3bc03254ca8997d0 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -14,7 +14,7 @@ Nrows = ceil(Nframes/4);
Ncols = ceil(Nframes/Nrows);
%
TNAME = [];
-FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
+FIGURE.fig = figure; %set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
for i_ = 1:numel(FRAMES)
frame_max = max(max(max(abs(toplot.FIELD(:,:,i_)))));
subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
diff --git a/matlab/plot/plot_param_evol.m b/matlab/plot/plot_param_evol.m
index e1f13b2cde404a860bc09bfd367e10b17d53db8b..6d16c6e37c8595169f189ef7585cfe69a26434aa 100644
--- a/matlab/plot/plot_param_evol.m
+++ b/matlab/plot/plot_param_evol.m
@@ -5,6 +5,7 @@ DT_EVOL = data.DT_EVOL;
CO_EVOL = data.CO_EVOL;
MUx_EVOL = data.MUx_EVOL;
MUy_EVOL = data.MUy_EVOL;
+MUz_EVOL = data.MUz_EVOL;
% K_T_EVOL = data.K_T_EVOL;
K_N_EVOL = data.K_N_EVOL;
L_EVOL = data.L_EVOL;
@@ -31,6 +32,7 @@ subplot(4,1,1);
subplot(4,1,2);
plot(TJOB_SE,MUx_EVOL,'DisplayName','$\mu_x$'); hold on;
plot(TJOB_SE,MUy_EVOL,'DisplayName','$\mu_y$'); hold on;
+ plot(TJOB_SE,MUz_EVOL,'DisplayName','$\mu_z$'); hold on;
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
xticks([]);
plot_tjob_lines(TJOB_SE,ylim)
diff --git a/matlab/setup.m b/matlab/setup.m
index 2c66ec4875c73062ca8b9a8340f58314a994c09e..468f08c57ea4feab7a6d2d3ce9319d126347fe6d 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -34,6 +34,7 @@ MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
MODEL.N_HD = N_HD;
MODEL.mu_z = MU_Z;
+MODEL.HYP_V = ['''',HYP_V,''''];
MODEL.mu_p = MU_P;
MODEL.mu_j = MU_J;
MODEL.nu = NU; % hyper diffusive coefficient nu for HW
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 92c9d1f4046e97c7b8673aeb73d77f3e38c6daf1..54e7f08b71adccef5c9afbd3bd596c64ed63fdb9 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -64,6 +64,7 @@ fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']);
fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']);
fprintf(fid,[' N_HD = ', num2str(MODEL.N_HD),'\n']);
fprintf(fid,[' mu_z = ', num2str(MODEL.mu_z),'\n']);
+fprintf(fid,[' HYP_V = ', MODEL.HYP_V,'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']);
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index fa1f1e5dd1518f45939de35b04cca87b6ec5a7e1..576e8ae1575819401fb4a35fc344e7ff19caebdd 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -14,7 +14,7 @@ implicit none
PRIVATE
! Geometry input parameters
CHARACTER(len=16), &
- PUBLIC, PROTECTED :: geom
+ PUBLIC, PROTECTED :: geom = 's-alpha'
REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
@@ -116,7 +116,7 @@ CONTAINS
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha geometry'
call eval_salpha_geometry
- CASE('Z-pinch')
+ CASE('Z-pinch','z-pinch','Zpinch','zpinch')
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
SHEARED = .FALSE.
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..e2108230318c033d3d6fadeb4e9950069dfab407
--- /dev/null
+++ b/testcases/smallest_problem/fort_00.90
@@ -0,0 +1,83 @@
+&BASIC
+ nrun = 10
+ dt = 0.01
+ tmax = 1
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 2
+ Lx = 200
+ Ny = 2
+ Ly = 60
+ Nz = 4
+ SG = .f.
+/
+&GEOMETRY
+ geom = 's-alpha'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ parallel_bc = 'dirichlet'
+/
+&OUTPUT_PAR
+ nsave_0d = 1
+ nsave_1d = -1
+ nsave_2d = -1
+ nsave_3d = 1
+ nsave_5d = 1
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .f.
+ write_Napj = .t.
+ write_Sapj = .f.
+ write_dens = .t.
+ write_temp = .t.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .t.
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ mu_p = 0
+ mu_j = 0
+ nu = 0.1
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ne = 2.0
+ K_Te = 0.4
+ K_Ni = 2.0
+ K_Ti = 0.4
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ gyrokin_CO = .true.
+ interspecies = .true.
+ !mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'phi'
+ ACT_ON_MODES = 'donothing'
+ init_background = 0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/smallest_problem/gyacomo_alpha b/testcases/smallest_problem/gyacomo_alpha
new file mode 120000
index 0000000000000000000000000000000000000000..0663323a521aaa271870a36110a656eb277a1f47
--- /dev/null
+++ b/testcases/smallest_problem/gyacomo_alpha
@@ -0,0 +1 @@
+../../bin/gyacomo_alpha
\ No newline at end of file
diff --git a/testcases/smallest_problem/gyacomo_debug b/testcases/smallest_problem/gyacomo_debug
new file mode 120000
index 0000000000000000000000000000000000000000..363e139a389f2e5d2d2097c09bf5ab363772dbfc
--- /dev/null
+++ b/testcases/smallest_problem/gyacomo_debug
@@ -0,0 +1 @@
+../../bin/gyacomo_debug
\ No newline at end of file
diff --git a/testcases/zpinch_example/gyacomo_debug b/testcases/zpinch_example/gyacomo_debug
new file mode 120000
index 0000000000000000000000000000000000000000..363e139a389f2e5d2d2097c09bf5ab363772dbfc
--- /dev/null
+++ b/testcases/zpinch_example/gyacomo_debug
@@ -0,0 +1 @@
+../../bin/gyacomo_debug
\ No newline at end of file
diff --git a/wk/CBC_P_J_scan.m b/wk/CBC_P_J_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..1b2a32eadf70e111f945ddb8ac49e1d852b43247
--- /dev/null
+++ b/wk/CBC_P_J_scan.m
@@ -0,0 +1,218 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_debug';
+EXECNAME = 'gyacomo';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 0;
+P_a = [2:20];
+J_a = [2:15];
+% P_a = [18:22];
+% J_a = [2:6];
+% collision setting
+CO = 'DG';
+HYP_V = 'dvpar4';
+NU = 0.02;
+K_T = 6.96;
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 5e-3;
+TMAX = 25;
+kymin = 0.3;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(J_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+j = 1;
+for P = P_a
+i = 1;
+for J = J_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % 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)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ DT = DT0/sqrt(J);
+ PMAXE = P; % Hermite basis size of electrons
+ JMAXE = J; % Laguerre "
+ PMAXI = P; % " ions
+ JMAXI = J; % "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ SPS0D = 1; % 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/2; % 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
+ 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
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; 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 = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
+ 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
+ if ~isempty(data_.NU_EVOL)
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+
+ [~,it1] = min(abs(data_.Ts3D-0.7*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 = data_.Ts3D(it1:it2);
+ 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_.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+ end
+j = j + 1;
+end
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+if 0
+%% Check time evolution
+figure;
+field = squeeze(sum(abs(data_.PHI),3)); field_t = squeeze(field(2,1,:)); to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+
+if 0
+%% Pcolor of the peak
+figure;
+[XX_,YY_] = meshgrid(J_a,P_a);
+% pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none'); axis ij;
+pclr=imagesc_custom(XX_,YY_,y_'.*(y_>0)');
+title(['$\nu_{',data.CO,'}=$',num2str(data.NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+xlabel('$\kappa_T$'); ylabel('$P$, $J=P/2$');
+colormap(bluewhitered)
+clb=colorbar;
+clb.Label.String = '$\gamma c_s/R$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
+%%
+if(numel(J_a)>1 && numel(P_a)>1)
+%% Save metadata
+ jmin = num2str(min(J_a)); jmax = num2str(max(J_a));
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_J_',jmin,'_',jmax,...
+ '_P_',pmin,'_',pmax,'_kT_',num2str(K_T),'_',CONAME,'_',num2str(NU),'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s1name = '$J$';
+metadata.s1 = J_a;
+metadata.s2name = '$P$';
+metadata.s2 = P_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+% metadata.input_file = h5read([LOCALDIR,'/outputs_00.h5'],'/files/STDIN.00');
+metadata.date = date;
+% tosave.data = metadata;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
\ No newline at end of file
diff --git a/wk/CBC_kT_scan_salpha.m b/wk/CBC_hypcoll_PJ_scan.m
similarity index 60%
rename from wk/CBC_kT_scan_salpha.m
rename to wk/CBC_hypcoll_PJ_scan.m
index 6229c864d14666e227f87eee981ea369af1f0778..1375c493d07f7b3a85d23e91cf269b5030f95e3a 100644
--- a/wk/CBC_kT_scan_salpha.m
+++ b/wk/CBC_hypcoll_PJ_scan.m
@@ -8,60 +8,64 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
-SIMID = 'p2_CBC_convergence_KT_PJ'; % Name of the simulation
-% SIMID = 'dbg'; % Name of the simulation
-RERUN = 0; % rerun if the data does not exist
-RUN = 0;
-KT_a = [3:0.5:7];
-% KT_a = [3];
-% P_a = [22];
+SIMID = 'p2_linear'; % Name of the simulation
+RERUN = 1; % If you want to rerun the sim (bypass the check of existing data)
+RUN = 1;
+SAVEDATA = 1;
+NU_a = [0.01:0.01:0.1 0.2 0.5 1.0];
P_a = [2:2:30];
-% P_a = [2:12];
+% NU_a = 0.8;
+% P_a = 8;
J_a = floor(P_a/2);
% collision setting
-CO = 'DG';
-NU = 0.05;
+% CO = 'none';
+CO = 'dvpar4';
+HYP_V = 'dvpar4';
+% HYP_V = 'hypcoll';
GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 1e-4; % electron plasma beta
+BETA = 0; % electron plasma beta
% background gradients setting
-K_N = 2.22; % Density '''
+K_N = 2.22;
+K_T = 6.96;
+% K_T = 5.3;
% Geometry
-% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT = 1e-3;
-TMAX = 25;
+DT0 = 1e-2;
+TMAX = 30;
kymin = 0.3;
NY = 2;
% arrays for the result
-g_ky = zeros(numel(KT_a),numel(P_a),NY/2+1);
+g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
j = 1;
for P = P_a
i = 1;
-for KT = KT_a
+for NU_ = NU_a
+ NU = NU_;
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
- K_Te = KT; % ele Temperature '''
- K_Ti = KT; % ion Temperature '''
- K_Ne = K_N; % ele Density '''
- K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
% P = 20;
J = floor(P/2);
+ DT = DT0/sqrt(J);
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
- NX = 12; % real space x-gridpoints
+ K_Ne = K_N; % ele Density '''
+ K_Te = K_T; % ele Temperature '''
+ K_Ni = K_N; % ion Density gradient drive
+ K_Ti = K_T; % ion Temperature '''
+ NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
@@ -82,7 +86,7 @@ for KT = KT_a
SPS0D = 1; % 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/2; % Sampling per time unit for 5D arrays
+ SPS5D = 1/10; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
@@ -96,11 +100,11 @@ for KT = KT_a
INIT_OPT= 'phi'; % 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_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;
+ W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -110,9 +114,9 @@ for KT = KT_a
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
- MU_Z = 0.2; %
- MU_P = 0.0; %
- MU_J = 0.0; %
+ MU_Z = 1.0; %
+ MU_P = NU; %
+ MU_J = NU; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
@@ -127,137 +131,121 @@ for KT = KT_a
data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
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'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
+ if ~isempty(data_.NU_EVOL)
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- % Load results after trying to run
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
-
- data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ [~,it1] = min(abs(data_.Ts3D-0.7*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 = data_.Ts3D(it1:it2);
+ 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,:));
- % 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 = 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;
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.ky(ikmax)); disp(msg);
+ end
end
- [gmax, ikmax] = max(g_ky(i,j,:));
-
- msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.ky(ikmax)); disp(msg);
-
-
+
i = i + 1;
-end
+ end
j = j + 1;
end
-if 0
+if 0
%% Check time evolution
figure;
+to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
-if 1
+%% take max growth rate among z coordinate
+[y_,idx_] = max(g_ky,[],3);
+e_ = g_std(:,:,idx_);
+
+if 1
%% Study of the peak growth rate
figure
-y_ = g_ky;
-e_ = 0.05;
-
-% filter growth rate with less than 0.05 value
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
-
-[y_,idx_] = max(g_ky,[],3);
-for i = 1:numel(idx_)
- e_ = g_std(:,:,idx_(i));
-end
-
-colors_ = jet(numel(KT_a));
+colors_ = jet(numel(NU_a));
subplot(121)
-for i = 1:numel(KT_a)
+for i = 1:numel(NU_a)
% errorbar(P_a,y_(i,:),e_(i,:),...
% 'LineWidth',1.2,...
-% 'DisplayName',['$\nu=$',num2str(KT_a(i))],...
+% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
% 'color',colors_(i,:));
plot(P_a,y_(i,:),'s-',...
'LineWidth',2.0,...
- 'DisplayName',['$\kappa_T=$',num2str(KT_a(i))],...
+ 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
'color',colors_(i,:));
hold on;
end
-title(['$\nu_{',CO,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=$',num2str(kymin)]);
legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
colors_ = jet(numel(P_a));
subplot(122)
for j = 1:numel(P_a)
-% errorbar(KT_a,y_(:,j),e_(:,j),...
-% 'LineWidth',1.2,...
-% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
-% 'color',colors_(j,:));
- plot(KT_a,y_(:,j),'s-',...
+ plot(NU_a,y_(:,j),'s-',...
'LineWidth',2.0,...
'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
'color',colors_(j,:));
hold on;
end
-title(['$\nu_{',CO,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
-legend('show'); xlabel('$\kappa_T$'); ylabel('$\gamma$');
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=$',num2str(kymin)]);
+legend('show'); xlabel(['$\nu_{',CO,'}$']); ylabel('$\gamma$');
end
if 0
%% Pcolor of the peak
figure;
-[XX_,YY_] = meshgrid(KT_a,P_a);
+[XX_,YY_] = meshgrid(NU_a,P_a);
% pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none'); axis ij;
pclr=imagesc_custom(XX_,YY_,y_'.*(y_>0)');
-title(['$\nu_{',data.CO,'}=$',num2str(data.NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
-xlabel('$\kappa_T$'); ylabel('$P$, $J=P/2$');
-colormap(bluewhitered)
+title(['$\kappa_T$',num2str(data_.K_T),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+xlabel('$\nu$'); ylabel('$P$, $J=P/2$');
+colormap(jet)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
-%%
+
+if(numel(NU_a)>1 && numel(P_a)>1)
%% Save metadata
-ktmin = num2str(min(KT_a)); ktmax = num2str(max(KT_a));
+numin = num2str(min(NU_a)); numax = num2str(max(NU_a));
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
- '_kT_',ktmin,'_',ktmax,'_',...
- '_P_',pmin,'_',pmax,'_',data_.CO,'_',num2str(NU),'.mat'];
+ '_nu_',numin,'_',numax,'_',HYP_V,...
+ '_P_',pmin,'_',pmax,'_KT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = kymin;
-metadata.title = ['$\nu_{',data_.CO,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
-metadata.par = data_.PARAMS;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
-metadata.s1name = '$\kappa_T$';
-metadata.s1 = KT_a;
-metadata.s2name = '$P$, $J=P/2$';
+metadata.s1name = ['$\nu_{',CONAME,'}$'];
+metadata.s1 = NU_a;
+metadata.s2name = '$P$, $J=\lfloor P/2 \rfloor$';
metadata.s2 = P_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
-metadata.input_file = h5read([data_.localdir,'/outputs_00.h5'],'/files/STDIN.00');
metadata.date = date;
-% tosave.data = metadata;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
-clear metadata tosave
\ No newline at end of file
+clear metadata tosave
+end
\ No newline at end of file
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 7ff5ff638efbed2e0050e71e250d69d4baab23ed..66912c279207326ed11e4447a6aefedb726a79d0 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -1,166 +1,278 @@
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-KT_a = [3:0.5:5];
-P_a_6 = [6 6 6 6 6 6 6];
-J_a_6 = [0 1 2 3 4 5 6];
-P_a_10 = 10*ones(1,6);
-J_a_10 = 5:10;
-% P_a_10b = 10*ones(1,7);
-% J_a_10b = 10:17;
-P_a_20 = 20*ones(1,11);
-J_a_20 = 10:20;
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_debug';
+EXECNAME = 'gyacomo';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+KT_a = [3:0.5:6.5 6.96];
+% KT_a = [5.5 6];
+% P_a = [25];
+P_a = [3:1:29];
+% P_a = [2:2:40];
+J_a = floor(P_a/2);
+% collision setting
+CO = 'DG';
+NU = 0.05;
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 1e-2;
+TMAX = 30;
+kymin = 0.3;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(KT_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+
+j = 1;
+for P = P_a
+i = 1;
+for KT = KT_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % 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)
+ K_Te = KT; % ele Temperature '''
+ K_Ti = KT; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+% P = 20;
+ J = floor(P/2);
+ DT = DT0/sqrt(J);
+ PMAXE = P; % Hermite basis size of electrons
+ JMAXE = J; % Laguerre "
+ PMAXI = P; % " ions
+ JMAXI = J; % "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % number of perpendicular planes (parallel grid)
+ NPOL = 1;
+ SG = 0; % Staggered z grids option
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ %% GEOMETRY
+ % GEOMETRY= 's-alpha';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ SPS0D = 1; % 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/2; % 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
+ 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
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% OUTPUTS
+ W_DOUBLE = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; 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 = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
+ 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
+ if ~isempty(data_.NU_EVOL)
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
-P_a = [P_a_20]; J_a = [J_a_20];
-% KT_a = 5.0; P_a = 20; J_a = 20;
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
-g_max= zeros(numel(P_a),numel(KT_a));
-g_avg= g_max*0;
-g_std= g_max*0;
-k_max= g_max*0;
+ % 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
-CO = 'DG'; GKCO = 0;
-NU = 0.0;
-DT = 7e-3;
-TMAX = 40;
-ky_ = 0.15;
-SIMID = 'linear_CBC_kT_threshold'; % Name of the simulation
-RUN = 0;
+ [~,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 = 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_.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
-for i = 1:numel(P_a)
-P = P_a(i); J = J_a(i);
- j=1;
- %Set Up parameters
- for K_Ti = KT_a
- CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
- TAU = 1.0; % e/i temperature ratio
- K_Ni = 2.22; K_Ne = K_Ni;
- K_Te = K_Ti; % Temperature '''
- SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
- KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
- BETA = 0e-1; % electron plasma beta
- PMAXE = P; JMAXE = J;
- PMAXI = P; JMAXI = J;
- NX = 8; % real space x-gridpoints
- NY = 6; % '' y-gridpoints
- LX = 2*pi/0.15; % Size of the squared frequency domain
- LY = 2*pi/ky_;
- NZ = 24; % number of perpendicular planes (parallel grid)
- NPOL = 1; SG = 0; NEXC = 1;
- GEOMETRY= 's-alpha';
- Q0 = 1.4; % safety factor
- SHEAR = 0.8; % magnetic shear (Not implemented yet)
- EPS = 0.18; % inverse aspect ratio
- SPS0D = 1; SPS2D = 0; SPS3D = 5;SPS5D= 1/5; SPSCP = 0;
- JOB2LOAD= -1;
- LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
- 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
- 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;
- 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 = 1.0; MU_P = 0.0; %
- MU_J = 0.0; LAMBDAD = 0.0;
- NOISE0 = 1.0e-4; % Init noise amplitude
- BCKGD0 = 0.0; % Init background
- k_gB = 1.0;k_cB = 1.0;
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
- %%-------------------------------------------------------------------------
- % RUN
- setup
- if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
- end
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
- % Load results
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [HELAZDIR,'results/',filename,'/'];
- data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+if 0
+%% Study of the peak growth rate
+figure
- %linear growth rate (adapted for 2D zpinch and fluxtube)
- trange = [0.5 1]*data.Ts3D(end);
- nplots = 0;
- 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);
+y_ = g_ky;
+e_ = 0.05;
- g_max(i,j) = gmax;
- k_max(i,j) = kmax;
-
- [g_avg(i,j), ik_] = max(lg.avg_g);
- g_std(i,j) = max(lg.std_g(ik_));
+% filter growth rate with less than 0.05 value
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
- j = j + 1;
- if 0
- %% Verify gamma time trace
- figure
- for ik_ = 1:numel(lg.ky)
- plot(lg.trange(2:end),lg.g_ky(ik_,:)','DisplayName',['$k_y=',num2str(lg.ky(ik_)),'$']); hold on;
- end
- xlabel('$t$'); ylabel('$\gamma$');
- title(data.param_title); legend('show');
- drawnow
- end
- end
+[y_,idx_] = max(g_ky,[],3);
+for i = 1:numel(idx_)
+ e_ = g_std(:,:,idx_(i));
end
-if 1
-%% PLOTS
-ERR_WEIGHT = 1/3; %weight of the error to compute marginal stability
-%% Superposed 1D plots
-sz_ = size(g_max);
-figure
-for i = 1:sz_(1)
- y_ = g_avg(i,:);
- e_ = g_std(i,:);
-
- y_ = y_.*(y_-e_*ERR_WEIGHT>0);
- e_ = e_ .* (y_>0);
- errorbar(KT_a,y_,e_,...
- 'LineWidth',1.2,...
- 'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']);
-% plot(KT_a,y_,...
+colors_ = jet(numel(KT_a));
+subplot(121)
+for i = 1:numel(KT_a)
+% errorbar(P_a,y_(i,:),e_(i,:),...
% 'LineWidth',1.2,...
-% 'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']);
+% 'DisplayName',['$\nu=$',num2str(KT_a(i))],...
+% 'color',colors_(i,:));
+ plot(P_a,y_(i,:),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['$\kappa_T=$',num2str(KT_a(i))],...
+ 'color',colors_(i,:));
hold on;
+end
+title(['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
+colors_ = jet(numel(P_a));
+subplot(122)
+for j = 1:numel(P_a)
+% errorbar(KT_a,y_(:,j),e_(:,j),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+% 'color',colors_(j,:));
+ plot(KT_a,y_(:,j),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+ 'color',colors_(j,:));
+ hold on;
+end
+title(['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+legend('show'); xlabel('$\kappa_T$'); ylabel('$\gamma$');
end
-title('Linear CBC $K_T$ threshold');
-legend('show'); xlabel('$K_T$'); ylabel('$\max_{k_y}(\gamma_k)$');
-drawnow
-%% Color map
-[NP__, KT__] = meshgrid(P_a+2*J_a, KT_a);
-% GG_ = g_avg;
-GG_ = g_avg .* (g_avg-g_std*ERR_WEIGHT > 0);
+
+if 0
+%% Pcolor of the peak
figure;
-% pclr = pcolor(KT__,NP__,g_max'); set(pclr,'EdgeColor','none');
-pclr = imagesc(KT_a,1:numel(P_a),GG_);
-LABELS = [];
-for i_ = 1:numel(P_a)
- LABELS = [LABELS; '(',sprintf('%2.0f',P_a(i_)),',',sprintf('%2.0f',J_a(i_)),')'];
+[XX_,YY_] = meshgrid(KT_a,P_a);
+% pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none'); axis ij;
+pclr=imagesc_custom(XX_,YY_,y_'.*(y_>0)');
+title(['$\nu_{',data.CO,'}=$',num2str(data.NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)]);
+xlabel('$\kappa_T$'); ylabel('$P$, $J=P/2$');
+colormap(bluewhitered)
+clb=colorbar;
+clb.Label.String = '$\gamma c_s/R$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
end
-yticks(1:numel(P_a));
-yticklabels(LABELS);
-xlabel('$\kappa_T$'); ylabel('$(P,J)$');
-title('Linear ITG threshold in CBC');
-colormap(bluewhitered);
%%
-%%
-end
-
+if(numel(KT_a)>1 && numel(P_a)>1)
+%% Save metadata
+ktmin = num2str(min(KT_a)); ktmax = num2str(max(KT_a));
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_kT_',ktmin,'_',ktmax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s1name = '$\kappa_T$';
+metadata.s1 = KT_a;
+metadata.s2name = '$P$, $J=\lfloor P/2 \rfloor$';
+metadata.s2 = P_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+% metadata.input_file = h5read([LOCALDIR,'/outputs_00.h5'],'/files/STDIN.00');
+metadata.date = date;
+% tosave.data = metadata;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
\ No newline at end of file
diff --git a/wk/CBC_nu_CO_scan_salpha.m b/wk/CBC_kT_nu_scan.m
similarity index 52%
rename from wk/CBC_nu_CO_scan_salpha.m
rename to wk/CBC_kT_nu_scan.m
index c4acfa75e14592313cc00f739eb2771298db2dcd..5bd21bfdbc4d2549113d5f148595982f848dc74c 100644
--- a/wk/CBC_nu_CO_scan_salpha.m
+++ b/wk/CBC_kT_nu_scan.m
@@ -5,59 +5,60 @@ addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_debug';
-% EXECNAME = 'gyacomo';
-EXECNAME = 'gyacomo_alpha';
+EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
-SIMID = 'convergence_CBC_salpha'; % Name of the simulation
-% SIMID = 'dbg'; % Name of the simulation
-RUN = 0; % to make sure it does not try to run
+SIMID = 'p2_linear'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
-
-% NU_a = [0.05];
-% P_a = [30];
-NU_a = [0.0 0.01 0.02 0.05 0.1];
-% P_a = [2 4:4:12];
-P_a = [2 4:4:28];
-J_a = floor(P_a/2);
+RUN = 1;
+KT_a = [3:0.5:6.5 6.96];
+NU_a = [0 0.01:0.01:0.1 0.2 0.5 1.0];
+% KT_a = 3.5;
+% NU_a = 0.5;
+P = 16;
+J = 8;
% collision setting
CO = 'DG';
GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
-KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
% background gradients setting
-K_Ne = 1*2.22; % ele Density '''
-K_Te = 1*6.96; % ele Temperature '''
-K_Ni = 1*2.22; % ion Density gradient drive
-K_Ti = 6.96; % ion Temperature '''
+K_N = 2.22; % Density '''
% Geometry
% GEOMETRY= 'miller';
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-% DT = 2e-3;
-DT = 5e-4;
-TMAX = 50;
-kymin = 0.4;
-NY = 2;
+DT = 2e-3;
+TMAX = 30;
+kymin = 0.3;
+NY = 2;
% arrays for the result
-g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
+g_ky = zeros(numel(KT_a),numel(NU_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
-
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+
j = 1;
-for P = P_a
-i = 1;
for NU = NU_a
+i = 1;
+for KT = KT_a
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% 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)
+ K_Te = KT; % ele Temperature '''
+ K_Ti = KT; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
%% GRID PARAMETERS
-% P = 20;
- J = floor(P/2);
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
@@ -97,11 +98,11 @@ for NU = NU_a
INIT_OPT= 'phi'; % 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_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;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -111,7 +112,7 @@ for NU = NU_a
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
- MU_Z = 0.2; %
+ MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
@@ -125,102 +126,84 @@ for NU = NU_a
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
- data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- if ((RERUN || isempty(data.NU_EVOL)) && RUN)
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
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
+ if ~isempty(data_.NU_EVOL)
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
- % Load results after trying to run
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
- data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- % 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
-% lg = compute_fluxtube_growth_rate(data,options);
-% [gmax, kmax] = max(lg.g_ky(:,end));
-% [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
-% g_ky(i,j,:) = g_ky;
-%
-% g_avg(i,j,:) = lg.avg_g;
-% g_std(i,j,:) = lg.std_g;
-
- [~,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 = 1: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);
- gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
- g_ky(i,j,ik) = gr(1);
- end
- [gmax, ikmax] = max(g_ky(i,j,:));
-
- msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
+ % 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 = 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_.ky(ikmax)); disp(msg);
+ end
+ end
i = i + 1;
end
j = j + 1;
end
-if 1
-%% Study of the peak growth rate
-figure
-
-y_ = g_ky;
-e_ = 0.05;
-
-% filter to noisy data
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
-
-[y_,idx_] = max(g_ky,[],3);
-for i = 1:numel(idx_)
- e_ = g_std(:,:,idx_(i));
-end
-
-colors_ = lines(numel(NU_a));
-subplot(121)
-for i = 1:numel(NU_a)
-% errorbar(P_a,y_(i,:),e_(i,:),...
-% 'LineWidth',1.2,...
-% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
-% 'color',colors_(i,:));
- plot(P_a,y_(i,:),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
- 'color',colors_(i,:));
- hold on;
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
-title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
-legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
-colors_ = jet(numel(P_a));
-subplot(122)
-for j = 1:numel(P_a)
-% errorbar(NU_a,y_(:,j),e_(:,j),...
-% 'LineWidth',1.2,...
-% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
-% 'color',colors_(j,:));
- plot(NU_a,y_(:,j),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
- 'color',colors_(j,:));
- hold on;
-end
-title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
-legend('show'); xlabel(['$\nu_{',CO,'}$']); ylabel('$\gamma$');
-end
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
-if 0
-%% Pcolor of the peak
-figure;
-[XX_,YY_] = meshgrid(NU_a,P_a);
-pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none');
+if(numel(KT_a)>1 && numel(NU_a)>1)
+%% Save metadata
+ktmin = num2str(min(KT_a)); ktmax = num2str(max(KT_a));
+numin = num2str(min(NU_a)); numax = num2str(max(NU_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
+ '_P_',num2str(P),'_J_',num2str(J),...
+ '_kT_',ktmin,'_',ktmax,...
+ '_nu_',numin,'_',numax,'_',CONAME,'.mat'];
+metadata.name = filename;
+metadata.kymin = kymin;
+metadata.title = ['P=',num2str(P),' J=',num2str(J),'$\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s1name = '$\kappa_T$';
+metadata.s1 = KT_a;
+metadata.s2name = ['$\nu_{',CONAME,'}$'];
+metadata.s2 = NU_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+metadata.date = date;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
end
\ No newline at end of file
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/CBC_ky_PJ_scan.m
deleted file mode 100644
index b53dde9de7c3fc399721803d27c6cf99d5ec2c47..0000000000000000000000000000000000000000
--- a/wk/CBC_ky_PJ_scan.m
+++ /dev/null
@@ -1,122 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-KY_a = 0.1:0.1:0.8;
-g_max= KY_a*0;
-g_avg= g_max*0;
-g_std= g_max*0;
-k_max= g_max*0;
-
-CO = 'DG'; GKCO = 0;
-NU = 0.05;
-DT = 2e-3;
-TMAX = 25;
-K_T = 6.96;
-SIMID = 'linear_CBC_circ_conv'; % Name of the simulation
-RUN = 0;
-% figure
-% P = 12;
-for P = [12]
-J = P/2;
-
-i=1;
-for ky_ = KY_a
-
-%Set Up parameters
-for j = 1
- CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
- TAU = 1.0; % e/i temperature ratio
- K_N = 2.22; K_Ne = K_N;
- K_Te = K_T; % Temperature '''
- SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
- KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
- BETA = 0e-1; % electron plasma beta
- PMAXE = P; JMAXE = J;
- PMAXI = P; JMAXI = J;
- NX = 12; % real space x-gridpoints
- NY = 2; % '' y-gridpoints
- LX = 2*pi/0.1; % Size of the squared frequency domain
- LY = 2*pi/ky_;
- NZ = 16; % number of perpendicular planes (parallel grid)
- NPOL = 1; SG = 0;
-% GEOMETRY= 's-alpha';
- GEOMETRY= 'circular';
- Q0 = 1.4; % safety factor
- SHEAR = 0.8; % magnetic shear (Not implemented yet)
- EPS = 0.18; % inverse aspect ratio
- SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
- JOB2LOAD= -1;
- LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
- 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
- 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;
- 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
-k_gB = 1.0;k_cB = 1.0;
-end
-%%-------------------------------------------------------------------------
-% RUN
-setup
-if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
-end
-
-% Load results
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
-data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
-
-%linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 0;
-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);
-
- g_max(i) = gmax;
- k_max(i) = kmax;
-
- g_avg(i) = lg.avg_g;
- g_std(i) = lg.std_g;
-
- i = i + 1;
-end
-%%
-
-% plot(KT_a,max(g_max,0));
-y_ = g_avg;
-e_ = g_std;
-
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
-% errorbar(KY_a,y_,e_,...
-% 'LineWidth',1.2,...
-% 'DisplayName',['(',num2str(P),',',num2str(J),')']);
-% hold on;
-% title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
-% legend('show'); xlabel('$k_y$'); ylabel('$\gamma$');
-% drawnow
-
-fig = plot_ballooning(data,options);
-
-end
-
diff --git a/wk/CBC_nu_CO_scan_miller.m b/wk/CBC_nu_CO_scan_miller.m
deleted file mode 100644
index e9a6384d6a67403907c816302eb206f29f60c639..0000000000000000000000000000000000000000
--- a/wk/CBC_nu_CO_scan_miller.m
+++ /dev/null
@@ -1,227 +0,0 @@
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % ... add
-addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
-addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
-addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
-% EXECNAME = 'gyacomo_dbg';
-% EXECNAME = 'gyacomo';
-EXECNAME = 'gyacomo_alpha';
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%
-% SIMID = 'linear_CBC_nu+PJ_scan_kT_6.96_SGGK'; % Name of the simulation
-SIMID = 'convergence_pITG_miller'; % Name of the simulation
-% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
-RUN = 0; % to make sure it does not try to run
-RERUN = 0; % rerun if the data does not exist
-
-% NU_a = [0.02];
-% P_a = [24];
-NU_a = [0.0 0.01 0.02 0.05 0.1];
-P_a = [2 4:4:28];
-% P_a = [24:4:28];
-J_a = floor(P_a/2);
-% collision setting
-CO = 'DG';
-GKCO = 0; % gyrokinetic operator
-COLL_KCUT = 1.75;
-% model
-KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 1e-4; % electron plasma beta
-% background gradients setting
-K_Ne = 1*2.22; % ele Density '''
-K_Te = 1*6.96; % ele Temperature '''
-K_Ni = 1*2.22; % ion Density gradient drive
-K_Ti = 6.96; % ion Temperature '''
-% Geometry
-GEOMETRY= 'miller';
-% GEOMETRY= 's-alpha';
-SHEAR = 0.8; % magnetic shear
-% time and numerical grid
-% DT = 2e-3;
-DT = 5e-4;
-TMAX = 50;
-kymin = 0.4;
-NY = 2;
-% arrays for the result
-g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
-g_avg= g_ky*0;
-g_std= g_ky*0;
-
-j = 1;
-for P = P_a
-i = 1;
-for NU = NU_a
- %% PHYSICAL PARAMETERS
- TAU = 1.0; % e/i temperature ratio
- % 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)
- %% GRID PARAMETERS
-% P = 20;
- J = floor(P/2);
- PMAXE = P; % Hermite basis size of electrons
- JMAXE = J; % Laguerre "
- PMAXI = P; % " ions
- JMAXI = J; % "
- NX = 8; % real space x-gridpoints
- LX = 2*pi/0.8; % Size of the squared frequency domain
- LY = 2*pi/kymin; % Size of the squared frequency domain
- NZ = 24; % number of perpendicular planes (parallel grid)
- NPOL = 1;
- SG = 0; % Staggered z grids option
- NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
- %% GEOMETRY
- % GEOMETRY= 's-alpha';
- EPS = 0.18; % inverse aspect ratio
- Q0 = 1.4; % safety factor
- KAPPA = 1.0; % elongation
- DELTA = 0.0; % triangularity
- ZETA = 0.0; % squareness
- % PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
- PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
- SHIFT_Y = 0.0;
- %% TIME PARMETERS
- SPS0D = 1; % 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/2; % 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
- 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
- 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.2; %
- MU_P = 0.0; %
- MU_J = 0.0; %
- LAMBDAD = 0.0;
- NOISE0 = 1.0e-5; % Init noise amplitude
- BCKGD0 = 0.0; % Init background
- k_gB = 1.0;
- k_cB = 1.0;
- %% RUN
- setup
- % naming
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
- % check if data exist to run if no data
- data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- if ((RERUN || isempty(data.NU_EVOL)) && RUN)
- 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 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
- end
-
- % Load results after trying to run
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
- data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
-
- % 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
-% lg = compute_fluxtube_growth_rate(data,options);
-% [gmax, kmax] = max(lg.g_ky(:,end));
-% [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
-% g_ky(i,j,:) = g_ky;
-%
-% g_avg(i,j,:) = lg.avg_g;
-% g_std(i,j,:) = lg.std_g;
-
- [~,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 = 1: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);
- gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
- g_ky(i,j,ik) = gr(1);
- end
- [gmax, ikmax] = max(g_ky(i,j,:));
-
- msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
-
-
- i = i + 1;
-end
-j = j + 1;
-end
-
-if 1
-%% Study of the peak growth rate
-figure
-
-y_ = g_ky;
-e_ = 0.05;
-
-% filter to noisy data
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
-
-[y_,idx_] = max(g_ky,[],3);
-for i = 1:numel(idx_)
- e_ = g_std(:,:,idx_(i));
-end
-
-colors_ = lines(numel(NU_a));
-subplot(121)
-for i = 1:numel(NU_a)
-% errorbar(P_a,y_(i,:),e_(i,:),...
-% 'LineWidth',1.2,...
-% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
-% 'color',colors_(i,:));
- plot(P_a,y_(i,:),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
- 'color',colors_(i,:));
- hold on;
-end
-title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
-legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
-
-colors_ = jet(numel(P_a));
-subplot(122)
-for j = 1:numel(P_a)
-% errorbar(NU_a,y_(:,j),e_(:,j),...
-% 'LineWidth',1.2,...
-% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
-% 'color',colors_(j,:));
- plot(NU_a,y_(:,j),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
- 'color',colors_(j,:));
- hold on;
-end
-title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
-legend('show'); xlabel(['$\nu_{',CO,'}$']); ylabel('$\gamma$');
-end
-
-if 0
-%% Pcolor of the peak
-figure;
-[XX_,YY_] = meshgrid(NU_a,P_a);
-pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none');
-end
\ No newline at end of file
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_PJ_scan.m
similarity index 73%
rename from wk/CBC_nu_CO_scan.m
rename to wk/CBC_nu_PJ_scan.m
index 437e0065fca4e069965fcf0cfa87d8a2deb9393d..8a28aca6ee1d1a177daf39a1a03329eddc6e1a4c 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_PJ_scan.m
@@ -8,14 +8,15 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
-SIMID = 'p2_CBC_convergence_coll_PJ'; % Name of the simulation
+SIMID = 'p2_linear'; % Name of the simulation
RERUN = 0; % If you want to rerun the sim (bypass the check of existing data)
RUN = 1;
-% NU_a = [0.0];
+% NU_a = [0.05];
% P_a = [8];
-NU_a = [0.0:0.01:0.1];
+NU_a = [0.01:0.01:0.1 0.2 0.5 1.0];
P_a = [2:2:30];
+% P_a = [32:2:40];
J_a = floor(P_a/2);
% collision setting
CO = 'DG';
@@ -23,17 +24,17 @@ GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 1e-4; % electron plasma beta
+BETA = 0; % electron plasma beta
% background gradients setting
K_N = 2.22;
-% K_T = 6.96;
K_T = 5.3;
+% K_T = 5.3;
% Geometry
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT = 1e-3;
-TMAX = 50;
+DT0 = 20e-3;
+TMAX = 60;
kymin = 0.3;
NY = 2;
% arrays for the result
@@ -52,6 +53,7 @@ for NU = NU_a
%% GRID PARAMETERS
% P = 20;
J = floor(P/2);
+ DT = DT0/sqrt(J);
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
@@ -60,7 +62,7 @@ for NU = NU_a
K_Te = K_T; % ele Temperature '''
K_Ni = K_N; % ion Density gradient drive
K_Ti = K_T; % ion Temperature '''
- NX = 12; % real space x-gridpoints
+ NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
@@ -74,14 +76,14 @@ for NU = NU_a
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
- % PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
- PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
SPS0D = 1; % 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/2; % Sampling per time unit for 5D arrays
+ SPS5D = 0; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
@@ -95,11 +97,11 @@ for NU = NU_a
INIT_OPT= 'phi'; % 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_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;
+ W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -109,7 +111,7 @@ for NU = NU_a
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
- MU_Z = 0.2; %
+ MU_Z = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
@@ -125,60 +127,55 @@ for NU = NU_a
% check if data exist to run if no data
data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
if RUN && (RERUN || isempty(data_.NU_EVOL) || numel(data_.Ts3D)<10)
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
+ if ~isempty(data_.NU_EVOL)
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
- % Load results after trying to run
- filename = [SIMID,'/',PARAMS,'/'];
- LOCALDIR = [gyacomodir,'results/',filename,'/'];
- data_ = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ [~,it1] = min(abs(data_.Ts3D-0.7*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 = data_.Ts3D(it1:it2);
+ 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,:));
- [~,it1] = min(abs(data_.Ts3D-0.8*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 = data_.Ts3D(it1:it2);
- 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;
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.ky(ikmax)); disp(msg);
+ end
end
- [gmax, ikmax] = max(g_ky(i,j,:));
-
- msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.ky(ikmax)); disp(msg);
-
-
+
i = i + 1;
-end
+ end
j = j + 1;
end
if 0
%% Check time evolution
figure;
-plot(data_.Ts3D,to_measure); hold on
-plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+plot(data_.Ts3D,log(field_t)); hold on
+plot(tw,to_measure,'--');
end
-if 1
-%% Study of the peak growth rate
-figure
-
-y_ = g_ky;
-e_ = 0.05;
-
-% filter to noisy data
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
+%% take max growth rate among z coordinate
[y_,idx_] = max(g_ky,[],3);
-for i = 1:numel(idx_)
- e_ = g_std(:,:,idx_(i));
-end
+e_ = g_std(:,:,idx_);
+
+if 1
+%% Study of the peak growth rate
+figure
colors_ = jet(numel(NU_a));
subplot(121)
@@ -228,23 +225,21 @@ end
numin = num2str(min(NU_a)); numax = num2str(max(NU_a));
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',num2str(kymin),...
- '_nu_',numin,'_',numax,'_',...
- '_P_',pmin,'_',pmax,'_KT_',num2str(K_Ti),'.mat'];
+ '_nu_',numin,'_',numax,'_',CONAME,...
+ '_P_',pmin,'_',pmax,'_KT_',num2str(K_T),'.mat'];
metadata.name = filename;
metadata.kymin = kymin;
-metadata.title = ['$\kappa_T$',num2str(K_Ti),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
-metadata.par = data_.PARAMS;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),', $\kappa_N=$',num2str(K_N),', $k_y=$',num2str(kymin)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
-metadata.s1name = '$\nu$';
+metadata.s1name = ['$\nu_{',CONAME,'}$'];
metadata.s1 = NU_a;
metadata.s2name = '$P$, $J=P/2$';
metadata.s2 = P_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
-metadata.input_file = h5read([data_.localdir,'/outputs_00.h5'],'/files/STDIN.00');
metadata.date = date;
-% tosave.data = metadata;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
-clear metadata tosave
+clear metadata tosave
\ No newline at end of file
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 09e9901ffed05d61cb4dfa88f8c67d146f393f3e..9ab388eac102be913b941cd9fa9b9594a3c8135c 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -57,6 +57,11 @@ folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
+% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
+
+% Miller CBC
+% folder = '/misc/gene_results/CBC/Miller_KT_6.96_128x64x24x16x8/';
+% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
% debug ? shearless
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_00/';
@@ -74,7 +79,7 @@ dashboard(gene_data);
%% Separated plot routines
if 0
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.1*gene_data.Ts3D(end);
+options.TAVG_0 = 0.5*gene_data.Ts3D(end);
options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
@@ -87,7 +92,7 @@ end
if 0
%% statistical transport averaging
-options.T = [100 500];
+options.T = [200 500];
fig = statistical_transport_averaging(gene_data,options);
end
@@ -109,7 +114,7 @@ options.NAME = '\phi';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'yz';
options.COMP = 'avg';
options.TIME = [50 200 500];
options.RESOLUTION = 256;
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index b05f7cb712da980254dbd53e62bfdcd566d9819b..a41a72ef2756e674895a8bcb8dd905c48a27ad5b 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -71,48 +71,48 @@ if 0
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}';
% options.NAME = 'n_i^{NZ}';
% options.NAME = 'Q_x';
-options.NAME = 'n_i';
+% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'yz';
+options.PLAN = 'kxz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [800:1000];
+options.TIME = [0:1500];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
create_film(data,options,'.gif')
end
-if 0
+if 1
%% fields snapshots
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
% options.NAME = 'T_i';
-% options.NAME = 'n_i-n_e';
+% options.NAME = 'n_i';
% options.NAME = '\phi^{NZ}';
-% options.NAME = 'N_i^{00}';
-% options.NAME = 'N_i^{00}-N_e^{00}';
+options.NAME = 'N_i^{00}';
+% options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
-options.COMP = 'avg';
-options.TIME = [50 200 500 1000];
+options.PLAN = 'kxky';
+options.COMP = 24;
+options.TIME = [600 700 800 1000];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
@@ -160,7 +160,7 @@ if 0
options.P2J = 0;
options.ST = 1;
options.NORMALIZED = 0;
-options.TIME = [180:10000];
+options.TIME = [500:800];
fig = show_moments_spectrum(data,options);
% fig = show_napjz(data,options);
% save_figure(data,fig,'.png');
@@ -202,7 +202,7 @@ end
if 0
%% Mode evolution
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
options.TIME = [000:9000];
options.KX_TW = [1 20]; %kx Growth rate time window
options.KY_TW = [0 20]; %ky Growth rate time window
diff --git a/wk/Ajay_scan_CH4_lin_ITG.m b/wk/benchmark scripts/Ajay_scan_CH4_lin_ITG.m
similarity index 100%
rename from wk/Ajay_scan_CH4_lin_ITG.m
rename to wk/benchmark scripts/Ajay_scan_CH4_lin_ITG.m
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 9a554ec4f1b801335fee2997acca98acfa8a52e3..32541c52ec039e8d379ba866695806eab689cf84 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -3,127 +3,71 @@ gyacomodir = '/home/ahoffman/gyacomo/';
% Partition of the computer where the data have to be searched
PARTITION = '/misc/gyacomo_outputs/';
% PARTITION = gyacomodir;
-%% 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/old/128x64x16x5x3';
-% resdir = 'CBC/96x96x16x3x2_Nexc_6';
-% resdir = 'CBC/128x96x16x3x2_Nexc_0';
-% resdir = 'CBC/192x96x24x13x7';
-
-% resdir = 'CBC/128x96x16x3x2_Nexc_0_periodic_chi';
-% resdir = 'CBC/64x32x16x3x2_Nexc_0_periodic_chi';
-% resdir = 'CBC/64x32x16x3x2_Nexc_0_Miller';
-% resdir = 'CBC/64x32x16x3x2_Nexc_1';
-% resdir = 'CBC/64x32x16x3x2_Nexc_2';
-% resdir = 'CBC/64x32x16x3x2_Nexc_3';
-% resdir = 'CBC/64x32x16x3x2_Nexc_3_change_dt';
-% resdir = 'CBC/64x32x16x3x2_Nexc_6';
-% resdir = 'CBC/64x32x16x3x2_Nexc_0';
-% resdir = 'CBC/64x32x16x3x2_Nexc_0_change_dt';
-% resdir = 'CBC/128x96x16x3x2_Nexc_0';
-% resdir = 'CBC/128x96x16x3x2_Nexc_0';
-
-% 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';
+%% CBC results with nuDGDK = 0.05
+% low resolution (Dirichlet)
+% resdir = 'paper_2_nonlinear/kT_6.96/3x2x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24'; %+ diff study
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24_dv4_diff';
+resdir = 'paper_2_nonlinear/kT_6.96/optimal_muz_nu_5x3x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x128x64x24';
+% low resolution (Cyclic)
+% resdir = 'paper_2_nonlinear/kT_6.96/3x2x128x64x24_cyclic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24_cyclic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x128x64x24_cyclic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x128x64x24_cyclic_BC';
+% low resolution (periodic)
+% resdir = 'paper_2_nonlinear/kT_6.96/3x2x128x64x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x128x64x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x128x64x24_periodic_BC';
+% high resolution (Dirichlet)
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x192x96x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x192x96x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x192x96x24';
+% high resolution (Cyclic)
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x192x96x24_cyclic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x192x96x24_cyclic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x192x96x24_cyclic_BC';
+% high resolution (periodic)
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x192x96x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x192x96x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/9x5x192x96x24_periodic_BC';
+% resdir = 'paper_2_nonlinear/kT_6.96/11x6x192x96x24_periodic_BC';
-% 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';
+% high z resolution (Dirichlet
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x32';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x192x96x64';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x96';
-% 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 = 'lin_CBC/128x64x16_5x3_Lx_7.854_Ly_125.6637_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_6.96_nu_5.0e-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/Salpha_GX_comparison';
-% resdir = 'Mandell_benchmark/5x3';
-% resdir = 'Mandell_benchmark/new_5x3';
-% resdir = 'Mandell_benchmark/new_5x3_more_diff';
-% resdir = 'Mandell_benchmark/7x5';
-% resdir = 'Mandell_benchmark/new_7x5';
-% resdir = 'Mandell_benchmark/new_7x5_more_diff';
-% % Paper 2 simulations
-% convergence CBC and Dimits regime
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_7x4x128x64x24';
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_3x2x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_5x3x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_7x4x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_9x5x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_6.96/CBC_21x11x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_11x6x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_3x2x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_5x3x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_7x4x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_9x5x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_11x6x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_13x7x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_15x8x128x64x24_Nexc_5';
-% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_17x9x128x64x24_Nexc_5';
-% Scan in kT
-% resdir = 'paper_2_nonlinear/kT_scan_DGGK_0.05/9x5x128x64x24';
-% resdir = 'paper_2_nonlinear/CBC_rerun/rerun_CBC_3x2x128x64x16';
-% resdir = 'paper_2_nonlinear/CBC_rerun/rerun_CBC_5x3x128x64x16';
-% resdir = 'dev/Napjz_spectrum';
-% resdir = 'dev/CBC_test';
+%% CBC results with nuDGDK = 0.01
+% resdir = 'paper_2_nonlinear/kT_6.96_nu_0.01/5x3x192x96x24';
-% Reruns
-% resdir = 'paper_2_nonlinear/kT_6.96_rerun/5x3x192x96x24';
-% resdir = 'paper_2_nonlinear/kT_6.96_rerun/7x4x192x96x24';
-% resdir = 'paper_2_nonlinear/kT_6.96_rerun/9x5x192x96x24';
-% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24';
-% resdir = 'paper_2_nonlinear/kT_6.96/7x4x128x64x24';
-% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/5x3x192x96x24';
-% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/7x4x192x96x24';
-% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/9x5x192x96x24';
+%% kT=5.3 results
+% resdir = 'paper_2_nonlinear/kT_5.3/5x3x128x64x64';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_MUxy_0';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_NL_-1';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_nuDG_0.01';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x192x96x64';
+% resdir = 'paper_2_nonlinear/kT_5.3/9x5x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x24';
-resdir = 'paper_2_nonlinear/kT_6.96/2jm1_version_5x3x128x64x24';
-% resdir = 'paper_2_nonlinear/dbg/CBC_test_Napjmir_3x2x128x64x24';
-% resdir = 'paper_2_nonlinear/dbg/CBC_bench_3x2x128x64x24';
+%% Old stuff
+% resdir = 'CBC/kT_4.5_128x64x16x13x7/';
+% resdir = 'CBC/kT_5.3_192x96x24x13x7/';
-% debug shearless
-% resdir = 'paper_2_nonlinear/shearless_CBC/7x4x192x96x24';
+%% Miller
+% resdir = 'paper_2_nonlinear/kT_4.15_miller/7x4x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_4.15_miller/7x4x128x64x32';
+% resdir = 'paper_2_nonlinear/kT_4.15_miller/7x4x128x64x64';
+% resdir = 'paper_2_nonlinear/kT_4.15_miller/7x5x192x64x24';
+% resdir = 'paper_2_nonlinear/kT_4.15_miller/GX_fig4_7x5x192x64x24';
-resdir = ['results/',resdir];
+%% Redo poster
+% resdir = 'paper_2_nonlinear/kT_6.96_Nexc1/5x3x128x64x24'; %% same as poster
+%%
JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_3D_Zpinch.m b/wk/lin_3D_Zpinch.m
index b56220c75ec9a62cec1906e096c8df8863027eca..e6ca8d59b90dca2d3c376562568d0485f8fdbfb1 100644
--- a/wk/lin_3D_Zpinch.m
+++ b/wk/lin_3D_Zpinch.m
@@ -171,12 +171,15 @@ save_figure(data,fig)
end
if 0
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
+options.NORMALIZED = 1;
+options.TIME = [000:9000];
+options.KX_TW = [1 20]; %kx Growth rate time window
+options.KY_TW = [0 20]; %ky Growth rate time window
options.NMA = 1;
options.NMODES = 1;
-options.iz = 'avg';
+options.iz = 'avg'; % avg or index
+options.ik = 1; % sum, max or index
+options.fftz.flag = 0;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
diff --git a/wk/lin_3Dzpinch.m b/wk/lin_3Dzpinch.m
deleted file mode 100644
index be6be1d9a9d3a4d7656d1b6cdbd7714107971a48..0000000000000000000000000000000000000000
--- a/wk/lin_3Dzpinch.m
+++ /dev/null
@@ -1,203 +0,0 @@
-%% 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"
-%
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % ... add
-addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
-addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
-addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
-SIMID = 'dbg'; % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-% EXECNAME = 'gyacomo_debug';
-% EXECNAME = 'gyacomo_alpha';
-EXECNAME = 'gyacomo';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.5; % Collision frequency
-TAU = 1e-2; % e/i temperature ratio
-K_Ne = 0; % ele Density '''
-K_Te = 0; % ele Temperature '''
-K_Ni = 0; % ion Density gradient drive
-K_Ti = 200; % ion Temperature '''
-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 = 2; % real space x-gridpoints
-NY = 40; % '' y-gridpoints
-LX = 2*pi/0.4; % Size of the squared frequency domain
-LY = 2*pi/0.2; % Size of the squared frequency domain
-NZ = 64; % number of perpendicular planes (parallel grid)
-NPOL = 2;
-SG = 0; % Staggered z grids option
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-%% GEOMETRY
-%% GEOMETRY
-GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-EPS = 0.18; % inverse aspect ratio
-Q0 = 1.4; % safety factor
-SHEAR = 0.0; % magnetic shear
-KAPPA = 1.0; % elongation
-DELTA = 0.0; % triangularity
-ZETA = 0.0; % squareness
-PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
-% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
-SHIFT_Y = 0.0;
-%% TIME PARMETERS
-TMAX = 20; % Maximal time unit
-DT = 1e-2; % Time step
-SPS0D = 1; % 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/2; % 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 = 'SG';
-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= 'phi'; % 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 = 1.0e-5; % Init noise amplitude
-BCKGD0 = 0.0; % Init background
-k_gB = 0.1;
-k_cB = 0.1;
-COLL_KCUT = 1.5;
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% Run linear simulation
-if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
-end
-
-%% Load results
-%%
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [gyacomodir,'results/',filename,'/'];
-FIGDIR = LOCALDIR;
-% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 01;
-data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
-
-%% Short analysis
-if 0
-%% linear growth rate (adapted for 2D zpinch and fluxtube)
-options.TRANGE = [0.5 1]*data.Ts3D(end);
-options.NPLOTS = 3; % 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 = [10 50];
-options.kymodes = [0.1 0.2 0.4];
-options.normalized = 1;
-% options.field = 'phi';
-fig = plot_ballooning(data,options);
-end
-
-if 0
-%% Hermite-Laguerre spectrum
-% options.TIME = 'avg';
-options.P2J = 0;
-options.ST = 0;
-options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';ls
-
-% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 1;
-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 1
-%% linear growth rate for 3D Zpinch (kz fourier transform)
-trange = [0.5 1]*data.Ts3D(end);
-options.INTERP = 1;
-options.kxky = 0;
-options.kzkx = 0;
-options.kzky = 1;
-[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
-end
-if 0
-%% ES Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
-% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
-options.NMA = 1;
-options.NMODES = 32;
-options.iz = 'avg';
-options.ik = 1;
-options.fftz.flag = 1;
-options.fftz.ky = 1.5; options.fftz.kx = 0;
-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
\ No newline at end of file
diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m
index fff23bffa8f194aba0001f96a3692323fb00a294..ed2d0e85e461944a284b79d2664b824ae187b232 100644
--- a/wk/lin_ETPY.m
+++ b/wk/lin_ETPY.m
@@ -13,33 +13,33 @@ SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_debug';
-EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
+EXECNAME = 'gyacomo_alpha';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
+NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 2.5; % ele Density '''
+K_Ne = 2.0; % 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)
+SIGMA_E = 1;%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 = 6;
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 = 40; % '' y-gridpoints
-LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.05; % Size of the squared frequency domain
+NX = 8; % real space x-gridpoints
+NY = 20; % '' y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain
+LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -47,10 +47,10 @@ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
%% GEOMETRY
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-EPS = 0.18; % inverse aspect ratio
-Q0 = 1.4; % safety factor
+EPS = 0.0; % inverse aspect ratio
+Q0 = 0.0; % safety factor
SHEAR = 0.0; % magnetic shear
-KAPPA = 1.0; % elongation
+KAPPA = 0.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
@@ -71,7 +71,7 @@ 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
@@ -80,11 +80,11 @@ KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % 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_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;
+W_DENS = 0; W_TEMP = 0;
+W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -94,7 +94,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 0.2; %
+MU_Z = 0.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
@@ -179,14 +179,15 @@ save_figure(data,fig)
end
if 1
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
-% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
+options.NORMALIZED = 1;
+options.TIME = [000:9000];
+options.KX_TW = [20 40]; %kx Growth rate time window
+options.KY_TW = [20 40]; %ky Growth rate time window
options.NMA = 1;
options.NMODES = 32;
-options.iz = 'avg';
-options.ik = 1;
+options.iz = 1; % avg or index
+options.ik = 1; % sum, max or index
+options.fftz.flag = 0;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 5390f57630162d2ef7dd64df3dd2a1a7d3e6c851..64bfc8880e7b1df9765dabcf53f1ee30eb8747fe 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -15,18 +15,17 @@ SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_debug';
-% EXECNAME = 'gyacomo_2jm1';
-% EXECNAME = 'gyacomo_dlnBdz';
-EXECNAME = 'gyacomo_alpha';
+EXECNAME = 'gyacomo';
+% EXECNAME = 'gyacomo_alpha';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.05; % Collision frequency
+NU = 0.0000; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 1*2.22; % ele Density '''
-K_Te = 1*6.96; % ele Temperature '''
+K_Ne = 0*2.22; % ele Density '''
+K_Te = 0*6.96; % ele Temperature '''
K_Ni = 1*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)
@@ -34,7 +33,7 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
+P = 6;
J = 2;%P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
@@ -49,21 +48,21 @@ NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
-% GEOMETRY= 's-alpha';
-GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
-PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
-% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+PARALLEL_BC = 'cyclic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
-TMAX = 40; % Maximal time unit
+TMAX = 25; % Maximal time unit
+% DT = 1e-2; % Time step
DT = 2e-2; % Time step
-% DT = 5e-4; % Time step
SPS0D = 1; % 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
@@ -98,8 +97,9 @@ 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_Z = 1.0; %
+HYP_V = 'dvpar4';
+MU_P = 0.5; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
@@ -113,8 +113,8 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'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,'/; time mpirun -np 4 ',gyacomodir,'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 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
diff --git a/wk/lin_RHT.m b/wk/lin_RHT.m
index 7c43f8fde0e8c47f26502d56ef0b230f4e8bfade..9fe1bcce57f6c9ebc80409997ecb5b3b35590902 100644
--- a/wk/lin_RHT.m
+++ b/wk/lin_RHT.m
@@ -15,34 +15,34 @@ EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.0001; % Collision frequency
TAU = 1.0; % e/i temperature ratio
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)
+% 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 = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.00; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
-J = P/2;
+P = 64;
+J = 16;
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 = pi/2.5; % Size of the squared frequency domain
+LX = 2*pi/0.05; % 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= 'miller';
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
@@ -52,11 +52,13 @@ NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) %0: adaptiv
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
-DT = 1e-2; % Time step
+% DT = 1e-2; % Time step
+DT = 1e-3; % Time step
SPS0D = 1; % 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
+SPS3D = 10; % Sampling per time unit for 2D arrays
+SPS5D = 1/10; % 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)
@@ -72,11 +74,11 @@ 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_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;
+W_DENS = 0; W_TEMP = 0;
+W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -101,26 +103,29 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% 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'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [GYACOMODIR,'results/',filename,'/'];
+LOCALDIR = [gyacomodir,'results/',filename,'/'];
+FIGDIR = LOCALDIR;
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
+JOBNUMMIN = 00; JOBNUMMAX = 01;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
%% 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);
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 3; % 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);
@@ -129,8 +134,8 @@ end
if 0
%% Ballooning plot
-options.time_2_plot = [120];
-options.kymodes = [0.1];
+options.time_2_plot = [10 50];
+options.kymodes = [0.1 0.2 0.4];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
@@ -142,9 +147,10 @@ if 0
options.P2J = 1;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-1D';ls
+
% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 0;
+options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [0 50];
options.specie = 'i';
@@ -157,33 +163,42 @@ end
if 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
+options.INTERP = 0;
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
-options.kzkx = 0;
-options.kzky = 0;
+options.kzkx = 1;
+options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
if 0
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 2;
-options.TIME = [0:50];
+options.NORMALIZED = 1;
+options.TIME = [000:9000];
+options.KX_TW = [20 40]; %kx Growth rate time window
+options.KY_TW = [20 40]; %ky Growth rate time window
options.NMA = 1;
-options.NMODES = 1;
-options.iz = 'avg';
+options.NMODES = 32;
+options.iz = 1; % avg or index
+options.ik = 1; % sum, max or index
+options.fftz.flag = 0;
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);
+ikx = 2; iky = 1; 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));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+t_ = data.Ts3D(it0:it1);
+TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.35*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2);
+clr_ = lines(20);
figure
-plot(data.Ts3D(it0:it1), plt(data.PHI));
+plot(t_, plt(data.PHI),'color',clr_(min((data.Pmaxi-1)/2,20),:)); hold on;
+plot(t_,0.5* exp(-t_*NU)+theory,'--k');
+plot([t_(1) t_(end)],theory*[1 1],'-k');
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
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index eeebd3e46b755b44ea99f4720fd43e3aa4cd152d..df8ed3ec4d00bb199394241ebf27e0f46f15cc6a 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -1,32 +1,51 @@
-%% Metadata locations
-% Scans over KT and PJ, keeping ky, CO constant
-% datafname = 'p2_CBC_convergence_KT_PJ/12x24_ky_0.3_kT_3_7__P_2_30_DGdkaa_0.05.mat';
-% datafname = 'p2_CBC_convergence_KT_PJ/12x24_ky_0.3_kT_3_7__P_2_30_DGdkaa_0.01.mat';
-% datafname = 'p2_CBC_convergence_KT_PJ/12x24_ky_0.3_kT_3_7__P_2_30_DGdkaa_0.mat';
-% Scans over NU and PJ, keeping ky and KY constant
-% datafname = 'p2_CBC_convergence_coll_PJ/12x24_ky_0.3_nu_0_0.1__P_2_30_KT_6.96.mat';
-datafname = 'p2_CBC_convergence_coll_PJ/12x24_ky_0.3_nu_0_0.1__P_2_30_KT_5.3.mat';
-
+% Metadata path
+%% Scans over KT and NU, keeping ky, CO constant (CBC_kT_nu_scan.m)
+% datafname = 'p2_linear/8x24_ky_0.3_P_4_J_2_kT_3_6.96_nu_0_1_DGDK.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_P_8_J_4_kT_3_6.96_nu_0_1_DGDK.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_P_16_J_8_kT_3_6.96_nu_0_1_DGDK.mat';
+%% Scans over KT and PJ, keeping ky, CO constant (CBC_kT_PJ_scan.m)
+% datafname = 'p2_linear/8x24_ky_0.3_kT_3_6.96_P_2_40_DGDK_0.05.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_kT_3_6.96_P_2_30_DGDK_0.05.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_kT_3_6.96_P_2_40_DGDK_0.025.mat';
+%% Scans over NU and PJ, keeping ky and KY constant (CBC_nu_PJ_scan.m)
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_DGDK_P_2_40_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_DGDK_P_2_40_KT_5.3.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_SGGK_P_2_20_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_hypcoll_P_2_30_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_hypcoll_P_2_30_KT_5.3.mat';
+ datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_30_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_10_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_6_10_KT_6.96.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_nu_0.01_1_dvpar4_P_2_10_KT_5.3.mat';
+%% Scans over P and J, keeping nu, ky and kT constant (CBC_P_J_scan.m)
+% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_30_kT_6.96_DGDK_0.05.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_30_kT_5.3_DGDK_0.05.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_20_kT_6.96_DGDK_0.02.mat';
+% datafname = 'p2_linear/8x24_ky_0.3_J_2_15_P_2_30_kT_5.3_DGDK_0.025.mat';
%% Load data
fname = ['../results/',datafname];
d = load(fname);
if 1
%% Pcolor of the peak
figure;
-[XX_,YY_] = meshgrid(d.s1,d.s2);
+% [XX_,YY_] = meshgrid(d.s1,d.s2);
+[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
+set(gca,'XTicklabel',d.s1)
+set(gca,'YTicklabel',d.s2)
+% colormap(jet)
colormap(bluewhitered)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
-if 0
+if 1
%% Scan along first dimension
figure
-colors_ = jet(numel(d.s2));
+colors_ = jet(numel(d.s2));
for i = 1:numel(d.s2)
% plot(d.s1,d.data(:,i),'s-',...
% 'LineWidth',2.0,...
@@ -42,7 +61,8 @@ xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
colormap(colors_);
clb = colorbar;
-caxis([d.s2(1),d.s2(end)]);
+caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
+clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
clb.YTick=d.s2;
clb.Label.String = d.s2name;
clb.Label.Interpreter = 'latex';
@@ -53,23 +73,25 @@ if 0
figure
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
-% plot(d.s2,d.data(i,:),'s-',...
-% 'LineWidth',2.0,...
-% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
-% 'color',colors_(i,:));
- errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
+ plot(d.s2,d.data(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
+% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
+% 'LineWidth',2.0,...
+% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
+% 'color',colors_(i,:));
hold on;
end
xlabel(d.s2name); ylabel(d.dname);title(d.title);
xlim([d.s2(1) d.s2(end)]);
-colormap(jet);
+colormap(jet(numel(d.s1)));
clb = colorbar;
-caxis([d.s1(1),d.s1(end)]);
+caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
+clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
+clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
@@ -77,24 +99,66 @@ end
if 0
%% Convergence analysis
figure
+% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
- target_ = d.data(i,end);
- eps_ = abs(target_ - d.data(i,1:end-1));
- semilogy(d.s2(1:end-1),eps_,'s',...
+% target_ = d.data(i,end);
+ target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
+ if target_ > 0
+ eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
+ semilogy(d.s2(1:end-1),eps_,'-s',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
hold on;
+ end
end
-xlabel(d.s2name); ylabel('$\epsilon$');title(d.title);
+xlabel(d.s2name); ylabel('$\epsilon_r$');title(d.title);
xlim([d.s2(1) d.s2(end)]);
-colormap(jet);
+colormap(colors_);
clb = colorbar;
-caxis([d.s1(1),d.s1(end)]);
+caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
+clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
+clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
grid on;
+end
+if 0
+%% Pcolor of the error
+figure; i_ = 0;
+% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
+% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
+target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
+% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
+% eps_ = log(abs(target_ - d.data)/abs(target_));
+eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
+sign_ = 1;%sign(d.data - target_);
+eps_ = d.data;
+for i = 1:numel(d.s1)
+ target_ = d.data(i,end);
+ eps_(i,:) = log(abs(target_ - d.data(i,1:end)));
+ if target_ > 0
+% eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
+% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
+% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
+ else
+ end
+end
+[XX_,YY_] = meshgrid(d.s1,d.s2);
+[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
+pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
+title(d.title);
+xlabel(d.s1name); ylabel(d.s2name);
+set(gca,'XTicklabel',d.s1)
+set(gca,'YTicklabel',d.s2)
+% colormap(jet)
+colormap(bluewhitered)
+% caxis([-10, 0]);
+clb=colorbar;
+clb.Label.String = '$\log(\epsilon_r)$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
end
\ No newline at end of file
diff --git a/wk/p2_heatflux_salpha_convergence.m b/wk/p2_heatflux_salpha_convergence.m
index 092f776fee7e526fbd9731428fbcfb75b644cf34..edeaa0afbc2289835c718dfe41619cfb55bdbf6d 100644
--- a/wk/p2_heatflux_salpha_convergence.m
+++ b/wk/p2_heatflux_salpha_convergence.m
@@ -1,18 +1,26 @@
%% Heat flux convergence for kt=6.96 and 5.3
figure
title('s-$\alpha$ turb. heat flux conv.');
-% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
-P = [2 4 12];
-Qx = [52.4 46.1 0];
-std_= [12.1 4.98 0];
+% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5 (cyclic BC)
+P = [2 4 6 8];
+Qx = [47.6 46.1 48.9 51.6];
+std_= [4.19 4.98 2.78 6.25];
errorbar(P,Qx,std_/2,'s-r',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
-% KT 6.96, nuDGDK = 0.05, 192x96x24, Nexc 5
+% % KT 6.96, nuDGDK = 0.05, 192x96x24, Nexc 5 (dirichlet BC)
+% P = [4 6 8 10 ];
+% Qx = [44.1 48.9 45.5 00.0];
+% std_= [9.63 6.13 13.8 0.00];
+% errorbar(P,Qx,std_/2,'o-r',...
+% 'LineWidth',2.0,'MarkerSize',8,...
+% 'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
+% xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
+% KT 6.96, nuDGDK = 0.05, 192x96x24, Nexc 5 (dirichlet BC)
P = [4 6 8 10 ];
-Qx = [42.9 44.9 00.0 00.0];
-std_= [9.15 7.12 0.00 0.00];
+Qx = [36.7 44.4 45.5 00.0];
+std_= [2.67 4.57 13.8 0.00];
errorbar(P,Qx,std_/2,'o-r',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 6.9, nuDGDK 0.05'); hold on