diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index b9b99f3ae976c17f8dc5a0735440b5f6c5b9667b..e64b370b5414d03f9f3bc4a9d21da581b81f0456 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -7,8 +7,8 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData];
end
-n0 = 1;
-n1 = 26;%numel(X_);
+n0 = 250;
+n1 = numel(X_);
figure;
mvm = @(x) movmean(x,1);
shift = X_(n0);
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 6f576a1dc7532418496bc130269965759b09d529..6940e1d1a07eb10a85eb821c2c7906803189268b 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -65,8 +65,11 @@ subplot(224)
%% FCGK
P_ = 4; J_ = 2;
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
-kp = 1.0;
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
+% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
+mat_file_name = '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5';
+
+kp = 2.0;
kp_a = h5read(mat_file_name,'/coordkperp');
[~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx);
@@ -93,7 +96,7 @@ subplot(224)
%% Single eigenvalue analysis
% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_4_J_2_N_50_kpm_4.0/scanfiles_00005/self.0.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_20_P_6_J_3_N_50_kpm_4.0/scanfiles_00042/self.0.h5';
+mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk.coulomb_NFLR_20_P_6_J_3_N_50_kpm_4.0/scanfiles_00042/self.0.h5';
matidx = 01;
@@ -116,12 +119,23 @@ if 0
%%
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
+ '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
+ '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG A', 'Hacked SG B'};
+CONAME_A = {'SG 20 10',...
+ 'PA 20 10',...
+ 'FC 10 5 NFLR 4',...
+ 'FC 10 5 NFLR 12',...
+ 'FC 10 5 NFLR 12 k<2', ...
+ 'FC 4 2 NFLR 6',...
+ 'FC 4 2 NFLR 12', ...
+ 'Hacked SG A',...
+ 'Hacked SG B'};
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
@@ -142,6 +156,7 @@ for j_ = 1:numel(mfns)
end
subplot(121)
legend('show'); grid on;
+ylim([0,100]);
xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)')
subplot(122)
legend('show'); grid on;
@@ -151,14 +166,14 @@ end
%% Van Kampen plot
if 0
%%
-kperp= 0.1;
+kperp= 1.5;
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/wk/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 10 5 NFLR 12'};
grow = {};
puls = {};
for j_ = 1:numel(mfns)
@@ -176,6 +191,6 @@ for j_ = 1:numel(mfns)
% plot(puls{j_}, grow{j_},'o','DisplayName',CONAME_A{j_}); hold on;
plot(grow{j_},'o','DisplayName',CONAME_A{j_}); hold on;
end
-legend('show'); grid on;
+legend('show'); grid on; title(['$k_\perp=$',num2str(kperp)]);
xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)')
end
\ No newline at end of file
diff --git a/matlab/setup.m b/matlab/setup.m
index db703d0145b228acf16d253ce3ae1f8acfedcc53..a11f45461138342c4220f61a4c568aa063465f70 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -62,7 +62,9 @@ switch CO
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD'
% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''';
- COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
+% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
+ COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
+% COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5''';
end
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 7352707ed92cc119d7498f66fa8ae6f541a307a0..2d4c4fffa44e3ae5ba0d869fd89d117a43eb317f 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -31,7 +31,7 @@ implicit none
! derivatives of magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
! Relative magnetic field strength
- REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL
! Relative strength of major radius
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
@@ -68,6 +68,9 @@ CONTAINS
call eval_1D_geometry
ELSE
SELECT CASE(geom)
+ CASE('circular')
+ IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
+ call eval_circular_geometry
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
@@ -143,7 +146,8 @@ CONTAINS
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
- hatB(iz,eo) = 1._dp / hatR(iz,eo)
+ hatB (iz,eo) = 1._dp / hatR(iz,eo)
+ hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
@@ -169,6 +173,71 @@ CONTAINS
!--------------------------------------------------------------------------------
!
+
+ SUBROUTINE eval_circular_geometry
+ ! evaluate circular geometry model
+ ! Ref: Lapilonne et al., PoP, 2009
+ implicit none
+ REAL(dp) :: X, kx, ky, Gamma1, Gamma2, Gamma3
+
+ parity: DO eo = 0,1
+ zloop: DO iz = izgs,izge
+ X = zarray(iz,eo) - eps*SIN(zarray(iz,eo)) ! chi = theta - eps sin(theta)
+
+ ! metric
+ gxx(iz,eo) = 1._dp
+ gxy(iz,eo) = shear*X - eps*SIN(X)
+ gxz(iz,eo) = - SIN(X)
+ gyy(iz,eo) = 1._dp + (shear*X)**2 - 2._dp*eps*COS(X) - 2._dp*shear*X*eps*SIN(X)
+ gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) - 2._dp*COS(X) - shear*X*SIN(X)
+ gzz(iz,eo) = 1._dp/eps**2 - 2._dp*COS(X)/eps
+ dxdR(iz,eo)= COS(X)
+ dxdZ(iz,eo)= SIN(X)
+
+ ! Relative strengh of radius
+ hatR(iz,eo) = 1._dp + eps*COS(X)
+ hatZ(iz,eo) = 1._dp + eps*SIN(X)
+
+ ! toroidal coordinates
+ Rc (iz,eo) = hatR(iz,eo)
+ phic(iz,eo) = X
+ Zc (iz,eo) = hatZ(iz,eo)
+
+ ! Jacobian
+ Jacobian(iz,eo) = q0*hatR(iz,eo)**2
+
+ ! Relative strengh of modulus of B
+ hatB (iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2)
+ hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
+
+ ! Derivative of the magnetic field strenght
+ gradxB(iz,eo) = -(COS(X) + eps*SIN(X)**2)/hatB(iz,eo)**2 ! Gene put a factor hatB^2 or 1/hatR^2 in this
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = eps * SIN(X) * (1._dp - eps*COS(X)) / hatB(iz,eo)**2 ! Gene put a factor hatB or 1/hatR in this
+
+ Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo)
+ Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)
+ Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)
+ ! Curvature operator
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky(iky, ikx, iz,eo) = Gamma1*((kx + shear*X*ky)*gradzB(iz,eo)/eps - gradxB(iz,eo)*ky) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
+
+ ENDDO zloop
+ ENDDO parity
+
+ END SUBROUTINE eval_circular_geometry
+ !
+ !--------------------------------------------------------------------------------
+ !
+
+
SUBROUTINE eval_zpinch_geometry
! evaluate s-alpha geometry model
implicit none
@@ -202,7 +271,8 @@ CONTAINS
Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
- hatB(iz,eo) = 1._dp
+ hatB (iz,eo) = 1._dp
+ hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
@@ -357,6 +427,7 @@ END SUBROUTINE set_ikx_zBC_map
CALL allocate_array( gradyB,izgs,izge, 0,1)
CALL allocate_array( gradzB,izgs,izge, 0,1)
CALL allocate_array( hatB,izgs,izge, 0,1)
+ CALL allocate_array( hatB_NL,izgs,izge, 0,1)
CALL allocate_array( hatR,izgs,izge, 0,1)
CALL allocate_array( hatZ,izgs,izge, 0,1)
CALL allocate_array( Rc,izgs,izge, 0,1)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 9d1626b254f8e40532cefce100422eb8f540b32c..2fe44a5454cc09fe7eace8128e92f65de54944cd 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -128,7 +128,7 @@ SUBROUTINE moments_eq_rhs_e
! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term
- - Sepj(ip,ij,iky,ikx,iz)
+ - hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
@@ -269,8 +269,8 @@ SUBROUTINE moments_eq_rhs_i
+ mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) &
! Collision term
+ TColl_i(ip,ij,iky,ikx,iz)&
- ! Nonlinear term
- - Sipj(ip,ij,iky,ikx,iz)
+ ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor
+ - hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..6d45978e37d91baf212cb7494426659e7c3572e1
--- /dev/null
+++ b/wk/CBC_kT_PJ_scan.m
@@ -0,0 +1,118 @@
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+KT_a = [3:1.0:7];
+g_max= KT_a*0;
+g_avg= KT_a*0;
+g_std= KT_a*0;
+k_max= KT_a*0;
+
+CO = 'DG'; GKCO = 0;
+NU = 0.01;
+DT = 4e-3;
+TMAX = 50;
+ky_ = 0.3;
+SIMID = 'linear_CBC_kT_scan_ky_0.3'; % Name of the simulation
+RUN = 0;
+figure
+P = 12;
+for J = [0 1 2 4 6]
+% J = P/2;
+
+i=1;
+for K_T = KT_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';
+ 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
+GRADB = 1.0;CURVB = 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(KT_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_T$'); ylabel('$\gamma$');
+drawnow
+end
+
diff --git a/wk/CBC_kT_PJ_scan_quick_run.m b/wk/CBC_kT_PJ_scan_quick_run.m
deleted file mode 100644
index 4245d7d5502f243fcbd62ea2bf150dfc5d65a54f..0000000000000000000000000000000000000000
--- a/wk/CBC_kT_PJ_scan_quick_run.m
+++ /dev/null
@@ -1,47 +0,0 @@
-KT_a = [3:0.5:7];
-g_max= KT_a*0;
-g_avg= KT_a*0;
-g_std= KT_a*0;
-k_max= KT_a*0;
-
-NU = 0.05;
-DT = 2e-3;
-TMAX = 50;
-ky_ = 0.3;
-SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation
-% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation
-RUN = 1;
-figure
-
-for P = [20]
-
-i=1;
-for K_T = KT_a
-
- quick_run
-
- 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(KT_a,y_,e_,...
- 'LineWidth',1.2,...
- 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
-hold on;
-title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
-legend('show'); xlabel('$K_T$'); ylabel('$\gamma$');
-drawnow
-end
-
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..949847b83d6628c21e9e32a3bfd42ca1a093e36a
--- /dev/null
+++ b/wk/CBC_nu_CO_scan.m
@@ -0,0 +1,116 @@
+% NU_a = [0.05 0.15 0.25 0.35 0.45];
+NU_a = [0:0.05:0.5];
+g_max= NU_a*0;
+g_avg= NU_a*0;
+g_std= NU_a*0;
+k_max= NU_a*0;
+CO = 'LD';
+
+K_T = 5.3;
+DT = 2e-3;
+TMAX = 30;
+ky_ = 0.3;
+SIMID = 'linear_CBC_nu_scan_kT_5.3_ky_0.3_LDGK'; % Name of the simulation
+RUN = 0;
+figure
+
+for P = [2 4 6]
+
+i=1;
+for NU = NU_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
+ J = P/2;
+ 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';
+ 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)
+ GKCO = 1; % gyrokinetic operator
+ ABCO = 1; % interspecies collisions
+ INIT_ZF = 0; ZF_AMP = 0.0;
+ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+ NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+ KERN = 0; % Kernel model (0 : GK)
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ 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
+GRADB = 1.0;CURVB = 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(NU_a,y_,e_,...
+ 'LineWidth',1.2,...
+ 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
+hold on;
+title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']);
+legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$');
+drawnow
+end
+
diff --git a/wk/CBC_nu_CO_scan_quick_run.m b/wk/CBC_nu_CO_scan_quick_run.m
deleted file mode 100644
index 6b9bc4fa9da0ccc5147004ef0518f97f5996685b..0000000000000000000000000000000000000000
--- a/wk/CBC_nu_CO_scan_quick_run.m
+++ /dev/null
@@ -1,48 +0,0 @@
-% NU_a = [0.05 0.15 0.25 0.35 0.45];
-NU_a = [0:0.05:0.5];
-g_max= NU_a*0;
-g_avg= NU_a*0;
-g_std= NU_a*0;
-k_max= NU_a*0;
-CO = 'LR';
-
-K_T = 5.3;
-DT = 2e-3;
-TMAX = 30;
-ky_ = 0.3;
-SIMID = 'linear_CBC_nu_scan_ky=0.3_CLOS_0_LRGK'; % Name of the simulation
-RUN = 1;
-figure
-
-for P = [2 4 6 10 12 20]
-
-i=1;
-for NU = NU_a
-
- quick_run
-
- 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(NU_a,y_,e_,...
- 'LineWidth',1.2,...
- 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
-hold on;
-title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']);
-legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$');
-drawnow
-end
-
diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m
index 96c8b4b25f0bab29082531886a86f23029bc0612..30103b346b8f1bd8ec959a0d19b3eee2de799060 100644
--- a/wk/Dimits_2000_fig3.m
+++ b/wk/Dimits_2000_fig3.m
@@ -1,27 +1,99 @@
-%192x96x16x3x2 simulation
-kT_Xi_GM_32 = ...
- [7.0 3.6;...
- 6.0 2.6;...
- 5.0 1.1;...
+%% Heat flux Qi [R/rhos^2/cs]
+kN = 2.22;
+%-------------- GM ---------------
+%192x96x16x3x2 kymin=0.05
+kT_Qi_GM_32 = ...
+ [...
+ 9.0 1.0e+2 4.3e+1;...
+ 7.0 3.6e+1 0;...
+ 6.0 2.6e+1 0;...
+ 5.0 1.1e+1 0;...
+ 4.5 7.8e+0 0;...
];
-%128x64x16x5x3 simulation
-kT_Xi_GM_53 = ...
- [7.0 4.6;...
- 5.3 1.9;...
- 5.0 1.5;...
+%128x64x16x5x3 kymin=0.05
+kT_Qi_GM_53 = ...
+ [...
+ 13. 1.3e+2 3.5e+1;...
+ 11. 9.7e+1 2.2e+1;...
+ 9.0 9.0e+1 2.8e+1;...
+ 7.0 4.6e+1 0;...
+ 5.3 1.9e+1 0;...
+ 5.0 1.5e+1 0;...
+ 4.5 9.7e+0 0;...9_128
];
-%128x64x16x24x12 simulation
-kT_Xi_GN = ...
- [7.0 5.0;...
- 5.3 2.4;...
- 4.5 nan;...
+%128x64x16x5x3 kymin=0.02 (large box)
+kT_Qi_GM_53_LB = ...
+ [...
+ 13. 1.1e+2 2.0e+1;...
];
-
+%-------------- GENE ---------------
+%128x64x16x24x12 kymin=0.05
+kT_Qi_GENE = ...
+ [...
+ 13. 2.0e+2 6.6e+1;...
+ 11. 3.3e+2 1.6e+2;...
+ 9.0 1.1e+2 0;...
+ 7.0 5.0e+1 0;...
+ 5.3 2.4e+1 0;...
+ 4.5 1.9e-1 0;...
+ ];
+%128x64x16x24x12 kymin=0.02 (large box)
+kT_Qi_GM_53_LB = ...
+ [...
+ 13. 2.7e+2 2.2e+1;...
+ 11. 1.9e+2 1.7e+1;...
+ ];
+%% Heat conductivity Xi [Ln/rhoi^2/cs] computed as Xi = Qi/kT/kN
+%init
+kT_Xi_GM_32 = kT_Qi_GM_32;
+kT_Xi_GM_53 = kT_Qi_GM_53;
+kT_Xi_GENE = kT_Qi_GENE;
+%scale
+for i = 2:3
+kT_Xi_GM_32(:,i) = kT_Qi_GM_32(:,i)./kT_Qi_GM_32(:,1)./kN;
+kT_Xi_GM_53(:,i) = kT_Qi_GM_53(:,i)./kT_Qi_GM_53(:,1)./kN;
+kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,1)./kN;
+end
+%% Dimits fig 3 data
+KT_DIM = [4.0 4.5 5.0 6.0 7.0 9.0 12. 14. 16. 18.];
+LLNL_GK_DIM = [5.0 0.0; ...
+ 7.0 2.5;...
+ 9.0 5.0;...
+ 12. 8.0;...
+ 14. 9.0;...
+ 16. 9.5;...
+ 18. 10.];
+UCOL_GK_DIM = [5.0 0.5;...
+ 6.0 1.0;...
+ 7.0 1.5];
+GFL__97_DIM = [4.0 4.0;...
+ 4.5 5.0;...
+ 5.0 6.5;...
+ 7.0 8.0;...
+ 9.0 11.];
+GFL__98_DIM = [5.0 1.5;...
+ 7.0 5.0;...
+ 9.0 7.5];
+%% Plot
figure;
-plot(kT_Xi_GM_32(:,1), kT_Xi_GM_32(:,2),'o-','DisplayName','192x96x16x3x2'); hold on
-plot(kT_Xi_GM_53(:,1), kT_Xi_GM_53(:,2),'o-','DisplayName','128x64x16x5x3'); hold on
-plot(kT_Xi_GN(:,1), kT_Xi_GN(:,2),'o--k','DisplayName','GENE 128x64x16x24x12');
-xlabel('$\kappa_T$'); ylabel('$\chi_i$');
+if 1
+xye = kT_Xi_GM_32;
+errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','192x96x16x3x2','LineWidth',2.0); hold on
+xye = kT_Xi_GM_53;
+errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','128x64x16x5x3','LineWidth',2.0); hold on
+xye = kT_Xi_GENE;
+errorbar(xye(:,1), xye(:,2),xye(:,3),'v-.','DisplayName','GENE 128x64x16x24x12','LineWidth',2.0);
+end
+if 1
+ plot(LLNL_GK_DIM(:,1),LLNL_GK_DIM(:,2),'dk--','DisplayName','Dimits GK, LLNL'); hold on
+ plot(UCOL_GK_DIM(:,1),UCOL_GK_DIM(:,2),'*k','DisplayName','Dimits PIC, U.COL');
+ plot(GFL__97_DIM(:,1),GFL__97_DIM(:,2),'+k','DisplayName','Dimits GFL 97');
+ plot(GFL__98_DIM(:,1),GFL__98_DIM(:,2),'sk','DisplayName','Dimits GFL 98');
+
+end
+plot([6.96 6.96],[0 12],'-.r','DisplayName','CBC');
+plot([4.0 4.00],[0 12],'-.b','DisplayName','$\kappa_T^{crit}$');
+xlabel('$\kappa_T$'); ylabel('$\chi_i$[$L_n/\rho_i^2 v_{thi}$]');
xlim([0 20]);
ylim([0 12]);
title('Dimits et al. 2000, Fig. 3');
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 1ba1f5aafbb5318c25dc6bba07ef9feb9d10322b..da93ddf3fbbc9ce60f40961ff0730d311901dbec 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -11,7 +11,6 @@ system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
data.FIGDIR = LOCALDIR;
@@ -23,7 +22,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2;
+options.TAVG_0 = 0.5*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2;
options.TAVG_1 = data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -50,13 +49,13 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xz';
+options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-% options.TIME = data.Ts3D(1:10:end);
-options.TIME = [200:2:400];
+options.TIME = data.Ts3D(1:2:end);
+% options.TIME = [550:2:750];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -65,7 +64,7 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
% options.NAME = '\phi';
@@ -74,11 +73,11 @@ options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [20 60 200 400 500];
+options.TIME = [550 650];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -139,7 +138,7 @@ options.NAME = 'N_i^{00}';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
-options.COMPXY = '2D'; % avg/sum/max/zero/ 2D plot otherwise
+options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 57ba009efb1863e599a7171de4b6abdd7680cd01..693e4723c7c53fafeb396108110145618ea4d0b5 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -12,10 +12,13 @@
% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
-folder = '/misc/gene_results/CBC/128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/196x96x20x32x16_00/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/';
+% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
+% folder = '/misc/gene_results/CBC/KT_13_128x64x16x24x12/';
+folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -41,7 +44,7 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.AXISEQUAL = 1;
+options.AXISEQUAL = 0;
% options.NAME = 'Q_x';
options.NAME = '\phi';
% options.NAME = 'n_i';
@@ -50,8 +53,8 @@ options.NAME = '\phi';
options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
-options.COMP = 1;
-options.TIME = [15 50 100 200];
+options.COMP = 'avg';
+options.TIME = [20 30 40 50 60];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -67,7 +70,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 45ec1e032f649ed9860c25e02d54ef67e9f1001e..4e2e073bc932024b432fc0d81d2e598ac260c4c7 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -43,8 +43,11 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'CBC/128x64x16x5x3';
% outfile = 'CBC/192x96x16x3x2';
% outfile = 'CBC/128x64x16x5x3';
+% outfile = 'CBC/kT_9_128x64x16x5x3';
+outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
+JOBNUMMIN = 00; JOBNUMMAX = 06;
-outfile = 'CBC/kT_scan_128x64x16x5x3';
+% outfile = 'CBC/kT_scan_128x64x16x5x3';
% outfile = 'CBC/kT_scan_192x96x16x3x2';
%% Linear CBC
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 8a7931f43108c194a9f2d02f1057d287d1b686cc..0769f58e818a82302ff60810bcc60bf650327214 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -3,7 +3,8 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
-% RUN = 1; % To run or just to load
+SIMID = 'Kinetic_ballooning_mode'; % Name of the simulation
+RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
@@ -14,40 +15,40 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-% NU = 0.0; % Collision frequency
+NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 2.22;%2.0; % ion Density gradient drive
-K_Ne = K_N; % ele Density gradient drive
-% K_T = 4.0;%0.25*K_N; % Temperature '''
+K_Ne = K_N; % ele Density gradient drive
+K_T = 6.96;%0.25*K_N; % Temperature '''
K_Te = K_T; % Temperature '''
-% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 0e-1; % electron plasma beta
+SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.03; % electron plasma beta
%% GRID PARAMETERS
-% P = 12;
+P = 8;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 12; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
+NY = 10; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-% LY = 2*pi/0.3; % Size of the squared frequency domain
-LY = 2*pi/ky_;
+LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
GEOMETRY= 's-alpha';
+% GEOMETRY= 'circular';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-% TMAX = 100; % Maximal time unit
-% DT = 2e-3; % Time step
+TMAX = 30; % Maximal time unit
+DT = 5e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
@@ -55,12 +56,10 @@ SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-% SIMID = 'linear_CBC'; % Name of the simulation
-% SIMID = 'scan'; % Name of the simulation
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';
+CO = 'DG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -78,7 +77,6 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
-kmax = NX*pi/LX;% Highest fourier mode
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
@@ -101,7 +99,7 @@ if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
@@ -117,7 +115,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
-nplots = 0;
+nplots = 1;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
diff --git a/wk/scan_quick_run.m b/wk/scan_quick_run.m
deleted file mode 100644
index 4245d7d5502f243fcbd62ea2bf150dfc5d65a54f..0000000000000000000000000000000000000000
--- a/wk/scan_quick_run.m
+++ /dev/null
@@ -1,47 +0,0 @@
-KT_a = [3:0.5:7];
-g_max= KT_a*0;
-g_avg= KT_a*0;
-g_std= KT_a*0;
-k_max= KT_a*0;
-
-NU = 0.05;
-DT = 2e-3;
-TMAX = 50;
-ky_ = 0.3;
-SIMID = 'linear_CBC_KT_scan_ky=0.3_CLOS_0_DGGK_0.05'; % Name of the simulation
-% SIMID = 'linear_CBC_PJ_scan_KT_4_ky=0.3_CLOS_0_TMAX_100'; % Name of the simulation
-RUN = 1;
-figure
-
-for P = [20]
-
-i=1;
-for K_T = KT_a
-
- quick_run
-
- 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(KT_a,y_,e_,...
- 'LineWidth',1.2,...
- 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
-hold on;
-title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
-legend('show'); xlabel('$K_T$'); ylabel('$\gamma$');
-drawnow
-end
-