Skip to content
Snippets Groups Projects
Commit f6cc280a authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

scripts updates

parent fef2ca96
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 469 additions and 214 deletions
matlab/extract_fig_data.m
+ 2
− 2
View file @ f6cc280a
...@@ -7,8 +7,8 @@ for i = 1:numel(dataObjs) ...@@ -7,8 +7,8 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData]; X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData]; Y_ = [Y_ dataObjs(i).YData];
end end
n0 = 1; n0 = 250;
n1 = 26;%numel(X_); n1 = numel(X_);
figure; figure;
mvm = @(x) movmean(x,1); mvm = @(x) movmean(x,1);
shift = X_(n0); shift = X_(n0);
......
matlab/plot/plot_cosol_mat.m
+ 23
− 8
View file @ f6cc280a
...@@ -65,8 +65,11 @@ subplot(224) ...@@ -65,8 +65,11 @@ subplot(224)
%% FCGK %% FCGK
P_ = 4; J_ = 2; P_ = 4; J_ = 2;
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/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
kp = 1.0; % 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'); kp_a = h5read(mat_file_name,'/coordkperp');
[~,matidx] = min(abs(kp_a-kp)); [~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx); matidx = sprintf('%5.5i',matidx);
...@@ -93,7 +96,7 @@ subplot(224) ...@@ -93,7 +96,7 @@ subplot(224)
%% Single eigenvalue analysis %% 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_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; matidx = 01;
...@@ -116,12 +119,23 @@ if 0 ...@@ -116,12 +119,23 @@ if 0
%% %%
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',... 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_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_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_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_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',... '/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 figure
for j_ = 1:numel(mfns) for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name); mat_file_name = mfns{j_}; disp(mat_file_name);
...@@ -142,6 +156,7 @@ for j_ = 1:numel(mfns) ...@@ -142,6 +156,7 @@ for j_ = 1:numel(mfns)
end end
subplot(121) subplot(121)
legend('show'); grid on; legend('show'); grid on;
ylim([0,100]);
xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)') xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)')
subplot(122) subplot(122)
legend('show'); grid on; legend('show'); grid on;
...@@ -151,14 +166,14 @@ end ...@@ -151,14 +166,14 @@ end
%% Van Kampen plot %% Van Kampen plot
if 0 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',... 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_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_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_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 = {}; grow = {};
puls = {}; puls = {};
for j_ = 1:numel(mfns) for j_ = 1:numel(mfns)
...@@ -176,6 +191,6 @@ 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(puls{j_}, grow{j_},'o','DisplayName',CONAME_A{j_}); hold on;
plot(grow{j_},'o','DisplayName',CONAME_A{j_}); hold on; plot(grow{j_},'o','DisplayName',CONAME_A{j_}); hold on;
end end
legend('show'); grid on; legend('show'); grid on; title(['$k_\perp=$',num2str(kperp)]);
xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)') xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)')
end end
\ No newline at end of file
matlab/setup.m
+ 3
− 1
View file @ f6cc280a
...@@ -62,7 +62,9 @@ switch CO ...@@ -62,7 +62,9 @@ switch CO
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''; COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD' 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_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 end
% Time integration and intialization parameters % Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4'''; TIME_INTEGRATION.numerical_scheme = '''RK4''';
......
src/geometry_mod.F90
+ 74
− 3
View file @ f6cc280a
...@@ -31,7 +31,7 @@ implicit none ...@@ -31,7 +31,7 @@ implicit none
! derivatives of magnetic field strength ! derivatives of magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
! Relative magnetic field strength ! Relative magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL
! Relative strength of major radius ! Relative strength of major radius
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients ! Some geometrical coefficients
...@@ -68,6 +68,9 @@ CONTAINS ...@@ -68,6 +68,9 @@ CONTAINS
call eval_1D_geometry call eval_1D_geometry
ELSE ELSE
SELECT CASE(geom) SELECT CASE(geom)
CASE('circular')
IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
call eval_circular_geometry
CASE('s-alpha') CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry' IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry call eval_salphaB_geometry
...@@ -143,7 +146,8 @@ CONTAINS ...@@ -143,7 +146,8 @@ CONTAINS
Jacobian(iz,eo) = q0*hatR(iz,eo) Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B ! 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 ! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
...@@ -169,6 +173,71 @@ CONTAINS ...@@ -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 SUBROUTINE eval_zpinch_geometry
! evaluate s-alpha geometry model ! evaluate s-alpha geometry model
implicit none implicit none
...@@ -202,7 +271,8 @@ CONTAINS ...@@ -202,7 +271,8 @@ CONTAINS
Jacobian(iz,eo) = 1._dp Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B ! 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 ! Derivative of the magnetic field strenght
gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this 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 ...@@ -357,6 +427,7 @@ END SUBROUTINE set_ikx_zBC_map
CALL allocate_array( gradyB,izgs,izge, 0,1) CALL allocate_array( gradyB,izgs,izge, 0,1)
CALL allocate_array( gradzB,izgs,izge, 0,1) CALL allocate_array( gradzB,izgs,izge, 0,1)
CALL allocate_array( hatB,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( hatR,izgs,izge, 0,1)
CALL allocate_array( hatZ,izgs,izge, 0,1) CALL allocate_array( hatZ,izgs,izge, 0,1)
CALL allocate_array( Rc,izgs,izge, 0,1) CALL allocate_array( Rc,izgs,izge, 0,1)
......
src/moments_eq_rhs_mod.F90
+ 3
− 3
View file @ f6cc280a
...@@ -128,7 +128,7 @@ SUBROUTINE moments_eq_rhs_e ...@@ -128,7 +128,7 @@ SUBROUTINE moments_eq_rhs_e
! Collision term ! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) & + TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term ! Nonlinear term
- Sepj(ip,ij,iky,ikx,iz) - hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) & IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
...@@ -269,8 +269,8 @@ SUBROUTINE moments_eq_rhs_i ...@@ -269,8 +269,8 @@ SUBROUTINE moments_eq_rhs_i
+ mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) & + mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) &
! Collision term ! Collision term
+ TColl_i(ip,ij,iky,ikx,iz)& + TColl_i(ip,ij,iky,ikx,iz)&
! Nonlinear term ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor
- Sipj(ip,ij,iky,ikx,iz) - hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) & IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
......
wk/CBC_kT_PJ_scan.m 0 → 100644
+ 118
− 0
View file @ f6cc280a
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
wk/CBC_kT_PJ_scan_quick_run.m deleted 100644 → 0
+ 0
− 47
View file @ fef2ca96
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
wk/CBC_nu_CO_scan.m 0 → 100644
+ 116
− 0
View file @ f6cc280a
% 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
wk/CBC_nu_CO_scan_quick_run.m deleted 100644 → 0
+ 0
− 48
View file @ fef2ca96
% 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
wk/Dimits_2000_fig3.m
+ 92
− 20
View file @ f6cc280a
%192x96x16x3x2 simulation %% Heat flux Qi [R/rhos^2/cs]
kT_Xi_GM_32 = ... kN = 2.22;
[7.0 3.6;... %-------------- GM ---------------
6.0 2.6;... %192x96x16x3x2 kymin=0.05
5.0 1.1;... 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 %128x64x16x5x3 kymin=0.05
kT_Xi_GM_53 = ... kT_Qi_GM_53 = ...
[7.0 4.6;... [...
5.3 1.9;... 13. 1.3e+2 3.5e+1;...
5.0 1.5;... 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 %128x64x16x5x3 kymin=0.02 (large box)
kT_Xi_GN = ... kT_Qi_GM_53_LB = ...
[7.0 5.0;... [...
5.3 2.4;... 13. 1.1e+2 2.0e+1;...
4.5 nan;...
]; ];
%-------------- 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; figure;
plot(kT_Xi_GM_32(:,1), kT_Xi_GM_32(:,2),'o-','DisplayName','192x96x16x3x2'); hold on if 1
plot(kT_Xi_GM_53(:,1), kT_Xi_GM_53(:,2),'o-','DisplayName','128x64x16x5x3'); hold on xye = kT_Xi_GM_32;
plot(kT_Xi_GN(:,1), kT_Xi_GN(:,2),'o--k','DisplayName','GENE 128x64x16x24x12'); errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','192x96x16x3x2','LineWidth',2.0); hold on
xlabel('$\kappa_T$'); ylabel('$\chi_i$'); 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]); xlim([0 20]);
ylim([0 12]); ylim([0 12]);
title('Dimits et al. 2000, Fig. 3'); title('Dimits et al. 2000, Fig. 3');
......
wk/analysis_HeLaZ.m
+ 8
− 9
View file @ f6cc280a
...@@ -11,7 +11,6 @@ system(['mkdir -p ',LOCALDIR]); ...@@ -11,7 +11,6 @@ system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD); system(CMD);
% Load outputs from jobnummin up to jobnummax % 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 = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR; data.localdir = LOCALDIR;
data.FIGDIR = LOCALDIR; data.FIGDIR = LOCALDIR;
...@@ -23,7 +22,7 @@ FMT = '.fig'; ...@@ -23,7 +22,7 @@ FMT = '.fig';
if 1 if 1
%% Space time diagramm (fig 11 Ivanov 2020) %% 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.TAVG_1 = data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
...@@ -50,13 +49,13 @@ options.NAME = '\phi'; ...@@ -50,13 +49,13 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}'; % options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
options.PLAN = 'xz'; options.PLAN = 'xy';
% options.NAME = 'f_i'; % options.NAME = 'f_i';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D(1:10:end); options.TIME = data.Ts3D(1:2:end);
options.TIME = [200:2:400]; % options.TIME = [550:2:750];
data.EPS = 0.1; data.EPS = 0.1;
data.a = data.EPS * 2000; data.a = data.EPS * 2000;
create_film(data,options,'.gif') create_film(data,options,'.gif')
...@@ -65,7 +64,7 @@ end ...@@ -65,7 +64,7 @@ end
if 0 if 0
%% 2D snapshots %% 2D snapshots
% Options % Options
options.INTERP = 0; options.INTERP = 1;
options.POLARPLOT = 0; options.POLARPLOT = 0;
options.AXISEQUAL = 0; options.AXISEQUAL = 0;
% options.NAME = '\phi'; % options.NAME = '\phi';
...@@ -74,11 +73,11 @@ options.NAME = 'N_i^{00}'; ...@@ -74,11 +73,11 @@ options.NAME = 'N_i^{00}';
% options.NAME = 'T_i'; % options.NAME = 'T_i';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e'; % options.NAME = 'k^2n_e';
options.PLAN = 'kxky'; options.PLAN = 'xy';
% options.NAME 'f_i'; % options.NAME 'f_i';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
options.TIME = [20 60 200 400 500]; options.TIME = [550 650];
data.a = data.EPS * 2e3; data.a = data.EPS * 2e3;
fig = photomaton(data,options); fig = photomaton(data,options);
% save_figure(data,fig) % save_figure(data,fig)
...@@ -139,7 +138,7 @@ options.NAME = 'N_i^{00}'; ...@@ -139,7 +138,7 @@ options.NAME = 'N_i^{00}';
options.PLAN = 'kxky'; options.PLAN = 'kxky';
options.COMPZ = 'avg'; options.COMPZ = 'avg';
options.OK = 0; 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.COMPT = 'avg';
options.PLOT = 'semilogy'; options.PLOT = 'semilogy';
fig = spectrum_1D(data,options); fig = spectrum_1D(data,options);
......
wk/analysis_gene.m
+ 8
− 5
View file @ f6cc280a
...@@ -12,10 +12,13 @@ ...@@ -12,10 +12,13 @@
% folder = '/misc/gene_results/LD_zpinch_1.6/'; % 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_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/'; % 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/196x96x20x32x16_00/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/'; % folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/'; % 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 = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data); gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1 if 1
...@@ -41,7 +44,7 @@ if 0 ...@@ -41,7 +44,7 @@ if 0
% Options % Options
options.INTERP = 1; options.INTERP = 1;
options.POLARPLOT = 0; options.POLARPLOT = 0;
options.AXISEQUAL = 1; options.AXISEQUAL = 0;
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
options.NAME = '\phi'; options.NAME = '\phi';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
...@@ -50,8 +53,8 @@ options.NAME = '\phi'; ...@@ -50,8 +53,8 @@ options.NAME = '\phi';
options.PLAN = 'xy'; options.PLAN = 'xy';
% options.NAME ='f_e'; % options.NAME ='f_e';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 1; options.COMP = 'avg';
options.TIME = [15 50 100 200]; options.TIME = [20 30 40 50 60];
gene_data.a = data.EPS * 2000; gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options); fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png') save_figure(gene_data,fig,'.png')
...@@ -67,7 +70,7 @@ options.NAME = '\phi'; ...@@ -67,7 +70,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}'; % options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x'; % options.NAME = '\Gamma_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
options.PLAN = 'xy'; options.PLAN = 'kxky';
% options.NAME = 'f_e'; % options.NAME = 'f_e';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
......
wk/header_3D_results.m
+ 4
− 1
View file @ f6cc280a
...@@ -43,8 +43,11 @@ helazdir = '/home/ahoffman/HeLaZ/'; ...@@ -43,8 +43,11 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'CBC/128x64x16x5x3'; % outfile = 'CBC/128x64x16x5x3';
% outfile = 'CBC/192x96x16x3x2'; % outfile = 'CBC/192x96x16x3x2';
% outfile = 'CBC/128x64x16x5x3'; % 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'; % outfile = 'CBC/kT_scan_192x96x16x3x2';
%% Linear CBC %% Linear CBC
......
wk/quick_run.m
+ 18
− 20
View file @ f6cc280a
...@@ -3,7 +3,8 @@ ...@@ -3,7 +3,8 @@
% from matlab framework. It is meant to run only small problems in linear % from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy" % 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 addpath(genpath('../matlab')) % ... add
default_plots_options default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/'; HELAZDIR = '/home/ahoffman/HeLaZ/';
...@@ -14,40 +15,40 @@ EXECNAME = 'helaz3'; ...@@ -14,40 +15,40 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS %% PHYSICAL PARAMETERS
% NU = 0.0; % Collision frequency NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio TAU = 1.0; % e/i temperature ratio
K_N = 2.22;%2.0; % ion Density gradient drive K_N = 2.22;%2.0; % ion Density gradient drive
K_Ne = K_N; % ele Density gradient drive K_Ne = K_N; % ele Density gradient drive
% K_T = 4.0;%0.25*K_N; % Temperature ''' K_T = 6.96;%0.25*K_N; % Temperature '''
K_Te = K_T; % Temperature ''' K_Te = K_T; % Temperature '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) 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) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0e-1; % electron plasma beta BETA = 0.03; % electron plasma beta
%% GRID PARAMETERS %% GRID PARAMETERS
% P = 12; P = 8;
J = P/2; J = P/2;
PMAXE = P; % Hermite basis size of electrons PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre " JMAXE = J; % Laguerre "
PMAXI = P; % " ions PMAXI = P; % " ions
JMAXI = J; % " JMAXI = J; % "
NX = 12; % real space x-gridpoints 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 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/0.1; % Size of the squared frequency domain
LY = 2*pi/ky_;
NZ = 16; % number of perpendicular planes (parallel grid) NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1; NPOL = 1;
SG = 0; % Staggered z grids option SG = 0; % Staggered z grids option
%% GEOMETRY %% GEOMETRY
% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps % GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
GEOMETRY= 's-alpha'; GEOMETRY= 's-alpha';
% GEOMETRY= 'circular';
Q0 = 1.4; % safety factor Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear (Not implemented yet) SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS %% TIME PARMETERS
% TMAX = 100; % Maximal time unit TMAX = 30; % Maximal time unit
% DT = 2e-3; % Time step DT = 5e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % 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 SPS3D = 1; % Sampling per time unit for 2D arrays
...@@ -55,12 +56,10 @@ SPS5D = 1/5; % Sampling per time unit for 5D arrays ...@@ -55,12 +56,10 @@ SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1; JOB2LOAD= -1;
%% OPTIONS %% OPTIONS
% SIMID = 'linear_CBC'; % Name of the simulation
% SIMID = 'scan'; % Name of the simulation
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator % Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
% CO = 'SG'; CO = 'DG';
GKCO = 1; % gyrokinetic operator GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0; INIT_ZF = 0; ZF_AMP = 0.0;
...@@ -78,7 +77,6 @@ W_NAPJ = 1; W_SAPJ = 0; ...@@ -78,7 +77,6 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused % unused
HD_CO = 0.0; % Hyper diffusivity cutoff ratio HD_CO = 0.0; % Hyper diffusivity cutoff ratio
kmax = NX*pi/LX;% Highest fourier mode
MU = 0.0; % Hyperdiffusivity coefficient MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; % MU_X = MU; %
...@@ -101,7 +99,7 @@ if RUN ...@@ -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,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',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 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']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end end
...@@ -117,7 +115,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from ...@@ -117,7 +115,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
if 1 if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube) %% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end); trange = [0.5 1]*data.Ts3D(end);
nplots = 0; nplots = 1;
lg = compute_fluxtube_growth_rate(data,trange,nplots); lg = compute_fluxtube_growth_rate(data,trange,nplots);
[gmax, kmax] = max(lg.g_ky(:,end)); [gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
......
wk/scan_quick_run.m deleted 100644 → 0
+ 0
− 47
View file @ fef2ca96
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment