diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index e410d05492312d249c2acb2acf5ef700eb280f0c..faa3c15219f218966060db36de1f713dba4cf3f4 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -64,11 +64,11 @@ subplot(224)
%% FCGK
-P_ = 6; J_ = 3;
-mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5';
-kp = 1.0; matidx = round(kp/(8/150));
-% matidx = sprintf('%5.5i',matidx);
-% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_1_kpm_4.0_NFLR_20_m_0.h5';
+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;
+kp_a = h5read(mat_file_name,'/coordkperp');
+[~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx);
FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
@@ -116,7 +116,8 @@ 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/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
};
CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 16'};
@@ -124,13 +125,14 @@ figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
kp_a = h5read(mat_file_name,'/coordkperp');
+% kp_a = kp_a(kp_a<=3);
gmax = zeros(size(kp_a));
wmax = zeros(size(kp_a));
for idx_ = 0:numel(kp_a)-1
matidx = sprintf('%5.5i',idx_);
MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
- gmax(idx_+1) = max(real(eig(MAT)));
- wmax(idx_+1) = max(imag(eig(MAT)));
+ gmax(idx_+1) = max(abs(real(eig(MAT))));
+ wmax(idx_+1) = max(abs(imag(eig(MAT))));
end
subplot(121)
plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
@@ -144,3 +146,35 @@ xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)')
legend('show'); grid on;
xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)')
end
+
+%% Van Kampen plot
+if 0
+%%
+kperp= 3.9;
+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/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
+ };
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 12 extended'};
+grow = {};
+puls = {};
+for j_ = 1:numel(mfns)
+ mat_file_name = mfns{j_}; disp(mat_file_name);
+ kp_a = h5read(mat_file_name,'/coordkperp');
+ [~,idx_] = min(abs(kp_a-kperp));
+ matidx = sprintf('%5.5i',idx_);
+ MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
+ grow{j_} = real(eig(MAT));
+ puls{j_} = imag(eig(MAT));
+end
+
+figure
+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;
+xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)')
+end
diff --git a/matlab/plot/plot_kernels.m b/matlab/plot/plot_kernels.m
index 400e63cc268c1ecbf3ef17cfeaac9e870ec396f3..a02b82cae6367eb22d4012c19859286eec4407cb 100644
--- a/matlab/plot/plot_kernels.m
+++ b/matlab/plot/plot_kernels.m
@@ -1,6 +1,6 @@
if 0
%% Kernels
-kmax=7;
+kmax=5;
nmax=6;
kx = linspace(0,kmax,100);
@@ -11,8 +11,9 @@ for n_ = 0:nmax
end
ylim_ = ylim;
plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
-plot(kx,J0,'-r','DisplayName','$J_0$');
-legend('show')
+% J0 = besselj(0,kx);
+% plot(kx,J0,'-r','DisplayName','$J_0$');
+legend('show'); xlabel('k'); title('Kernel functions $\mathcal{K}_n$');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -52,7 +53,7 @@ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% from Sum_n Kn x Ln x Lj to Sum_n Kn Sum_s dnjs x Ls
-vperp = [5];
+vperp = [2];
kx = linspace(0,10,200);
Jmax = 15;
j = 10;
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 9bdcd861187e940466122d53b61650a74df73bfb..7565d5fb3d4dab46a806aecf12189aece0d03ede 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -24,7 +24,7 @@ PUBLIC :: collision_readinputs, coll_outputinputs
PUBLIC :: compute_TColl
PUBLIC :: compute_lenard_bernstein, compute_dougherty
PUBLIC :: LenardBernstein_e, LenardBernstein_i!, LenardBernstein GK
-PUBLIC :: DoughertyGK_e, DoughertyGK_i!, Dougherty GK
+PUBLIC :: DoughertyGK_ee, DoughertyGK_ii!, Dougherty GK
PUBLIC :: load_COSOlver_mat, compute_cosolver_coll
PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i
@@ -173,26 +173,80 @@ CONTAINS
IF (KIN_E) THEN
DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e
DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
- CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl_)
- TColl_e(ip,ij,ikx,iky,iz) = TColl_
+ IF(gyrokin_CO) THEN
+ CALL DoughertyGK_ee(ip,ij,ikx,iky,iz,TColl_)
+ ELSE
+ CALL DoughertyDK_ee(ip,ij,ikx,iky,iz,TColl_)
+ ENDIF
+ TColl_e(ip,ij,ikx,iky,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
ENDIF
DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i
DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
- CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl_)
- TColl_i(ip,ij,ikx,iky,iz) = TColl_
+ IF(gyrokin_CO) THEN
+ CALL DoughertyGK_ii(ip,ij,ikx,iky,iz,TColl_)
+ ELSE
+ CALL DoughertyDK_ii(ip,ij,ikx,iky,iz,TColl_)
+ ENDIF
+ TColl_i(ip,ij,ikx,iky,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
END SUBROUTINE compute_dougherty
-
- SUBROUTINE DoughertyGK_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
- USE fields, ONLY: moments_e, phi
- USE grid, ONLY: parray_e, jarray_e, kxarray, kyarray, kparray, Jmaxe, ip0_e, ip1_e, ip2_e
- USE array, ONLY: kernel_e
- USE basic
- USE model, ONLY: sigmae2_taue_o2, qe_taue, nu_ee, sqrt_sigmae2_taue_o2
- USE time_integration, ONLY : updatetlevel
+ !******************************************************************************!
+ !! Doughtery driftkinetic collision operator for electrons
+ !******************************************************************************!
+ SUBROUTINE DoughertyDK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ COMPLEX(dp), INTENT(OUT) :: TColl_
+ COMPLEX(dp) :: upar, Tpar, Tperp
+ REAL(dp) :: j_dp, p_dp
+ !** Auxiliary variables **
+ p_dp = REAL(parray_e(ip_),dp)
+ j_dp = REAL(jarray_e(ij_),dp)
+ !** Assembling collison operator **
+ TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
+ IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10
+ TColl_ = TColl_ + nadiab_moments_e(ip1_e,1,ikx_,iky_,iz_)
+ ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20
+ TColl_ = TColl_ + twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_) &
+ - SQRT2*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_)
+ ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01
+ TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_) &
+ - SQRT2*twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_)
+ ENDIF
+ TColl_ = nu_ee * TColl_
+ END SUBROUTINE DoughertyDK_ee
+ !******************************************************************************!
+ !! Doughtery driftkinetic collision operator for ions
+ !******************************************************************************!
+ SUBROUTINE DoughertyDK_ii(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ COMPLEX(dp), INTENT(OUT) :: TColl_
+ COMPLEX(dp) :: upar, Tpar, Tperp
+ REAL(dp) :: j_dp, p_dp
+ !** Auxiliary variables **
+ p_dp = REAL(parray_i(ip_),dp)
+ j_dp = REAL(jarray_i(ij_),dp)
+ !** Assembling collison operator **
+ TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
+ IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p1j0
+ TColl_ = TColl_ + nadiab_moments_i(ip1_i,1,ikx_,iky_,iz_)
+ ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p2j0
+ TColl_ = TColl_ + twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_) &
+ - SQRT2*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_)
+ ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! kronecker p0j1
+ TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_) &
+ - SQRT2*twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_)
+ ENDIF
+ TColl_ = nu_i * TColl_
+ END SUBROUTINE DoughertyDK_ii
+ !******************************************************************************!
+ !! Doughtery gyrokinetic collision operator for electrons
+ !******************************************************************************!
+ SUBROUTINE DoughertyGK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -214,12 +268,10 @@ CONTAINS
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
! -nuee (p + 2j + b^2/2) Nepj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*moments_e(ip_,ij_,ikx_,iky_,iz_,updatetlevel)
+ TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
- ! Get adiabatic moment
- TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_,iz_,eo_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -231,13 +283,13 @@ CONTAINS
Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_e(ip0_e,in_ ,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Perpendicular velocity
uperp_ = uperp_ + be_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_e(ip2_e,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
@@ -255,7 +307,7 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxe+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_,eo_) * moments_e(ip1_e,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_,eo_) * nadiab_moments_e(ip1_e,in_,ikx_,iky_,iz_)
ENDDO
TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -273,11 +325,11 @@ CONTAINS
Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_e(ip0_e,in_,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_e(ip2_e,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
@@ -289,11 +341,11 @@ CONTAINS
! Multiply by electron-electron collision coefficient
TColl_ = nu_ee * TColl_
- END SUBROUTINE DoughertyGK_e
+ END SUBROUTINE DoughertyGK_ee
!******************************************************************************!
!! Doughtery gyrokinetic collision operator for ions
!******************************************************************************!
- SUBROUTINE DoughertyGK_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE DoughertyGK_ii(ip_,ij_,ikx_,iky_,iz_,TColl_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -316,12 +368,10 @@ CONTAINS
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
! -nui (p + 2j + b^2/2) Nipj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*moments_i(ip_,ij_,ikx_,iky_,iz_,updatetlevel)
+ TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
- ! Get adiabatic moment
- TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_,iz_,eo_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -333,17 +383,17 @@ CONTAINS
Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_i(ip0_i,in_ ,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Perpendicular velocity
uperp_ = uperp_ + bi_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_i(ip2_i,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
! Add energy restoring term
TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikx_,iky_,iz_,eo_)
TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_,eo_)
@@ -357,7 +407,7 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxi+1
! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_,eo_) * moments_i(ip1_i,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_,eo_) * nadiab_moments_i(ip1_i,in_,ikx_,iky_,iz_)
ENDDO
TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -375,15 +425,15 @@ CONTAINS
Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_i(ip0_i,in_,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_i(ip2_i,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -391,7 +441,7 @@ CONTAINS
! Multiply by ion-ion collision coefficient
TColl_ = nu_i * TColl_
- END SUBROUTINE DoughertyGK_i
+ END SUBROUTINE DoughertyGK_ii
!******************************************************************************!
!! compute the collision terms in a (Np x Nj x Nkx x Nky) matrix all at once
diff --git a/src/inital.F90 b/src/inital.F90
index c103b9149fe3a8c5eb01dec4e02db36082924224..e17b69f455a414ab0136fae8dd42a7b7e54114e8 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -24,6 +24,7 @@ SUBROUTINE inital
IF ( job2load .GE. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
CALL load_moments ! get N_0
+
CALL poisson ! compute phi_0=phi(N_0)
! through initialization
ELSE
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 797267b88d2c35d6f70f3d551b4eb40fc1125dfb..cc69a2b1845375eb87bffede814f6669edf55335 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -91,6 +91,7 @@ CONTAINS
!******************************************************************************!
SUBROUTINE load_output_adapt_pj
IMPLICIT NONE
+ INTEGER :: pmaxloop_e, pmaxloop_i, jmaxloop_e, jmaxloop_i, Nkx_cp, Nky_cp, Nz_cp
! Checkpoint filename
WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',job2load,'.h5'
@@ -98,6 +99,10 @@ CONTAINS
IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Resume from ", rstfile
! Open file
CALL openf(rstfile, fidrst,mpicomm=comm0)
+ ! Get grid info
+ CALL getatt(fidrst,"/data/input/" , "Nkx", Nkx_cp)
+ CALL getatt(fidrst,"/data/input/" , "Nky", Nky_cp)
+ CALL getatt(fidrst,"/data/input/" , "Nz", Nz_cp)
!!!!!!!!! Load electron moments
IF (KIN_E) THEN
! Get the checkpoint moments degrees to allocate memory
@@ -107,7 +112,8 @@ CONTAINS
IF (my_id .EQ. 0) WRITE(*,*) "Je_cp = ", jmaxe_cp
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
! Allocate the required size to load checkpoints moments
- CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
+ ! CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, 1,Nkx_cp, 1,Nky_cp, 1,Nz_cp)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! start with moments_e/000001
@@ -118,21 +124,28 @@ CONTAINS
n_ = n_ - 1 ! n_ is not a file so take the previous one n_-1
! Read state of system from checkpoint file and load every moment to change the distribution
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
- CALL getarrnd(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
+ ! CALL getarrnd(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
+ CALL getarr(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, 1:Nkx_cp, 1:Nky_cp, 1:Nz_cp))
! Initialize simulation moments array with checkpoints ones
! (they may have a larger number of polynomials, set to 0 at the begining)
moments_e = 0._dp;
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz = izs,ize
- moments_e(ip,ij,ikx,iky,iz,:) = moments_e_cp(ip,ij,ikx,iky,iz)
+ pmaxloop_e = min(ipe_e,pmaxe_cp+1)
+ jmaxloop_e = min(ije_e,jmaxe_cp+1)
+ IF (ips_e .LE. pmaxe_cp+1) THEN
+ DO ip=ips_e,pmaxloop_e
+ IF (ijs_e .LE. jmaxe_cp+1) THEN
+ DO ij=ijs_e,jmaxloop_e
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz = izs,ize
+ moments_e(ip,ij,ikx,iky,iz,:) = moments_e_cp(ip,ij,ikx,iky,iz)
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
+ ENDIF
+ ENDDO
+ ENDIF
! Deallocate checkpoint arrays
DEALLOCATE(moments_e_cp)
ENDIF
@@ -144,7 +157,8 @@ CONTAINS
IF (my_id .EQ. 0) WRITE(*,*) "Ji_cp = ", jmaxi_cp
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
! Allocate the required size to load checkpoints moments
- CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
+ ! CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, 1,Nkx_cp, 1,Nky_cp, 1,Nz_cp)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! start with moments_e/000001
@@ -156,21 +170,29 @@ CONTAINS
! Read state of system from checkpoint file and load every moment to change the distribution
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
- CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
+ ! CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
+ CALL getarr(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, 1:Nkx_cp, 1:Nky_cp, 1:Nz_cp))
+
! Initialize simulation moments array with checkpoints ones
! (they may have a larger number of polynomials, set to 0 at the begining)
- moments_i = 0._dp
- DO ip=1,pmaxi_cp+1
- DO ij=1,jmaxi_cp+1
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz = izs,ize
- moments_i(ip,ij,ikx,iky,iz,:) = moments_i_cp(ip,ij,ikx,iky,iz)
+ moments_i = 0._dp;
+ pmaxloop_i = min(ipe_i,pmaxi_cp+1)
+ jmaxloop_i = min(ije_i,jmaxi_cp+1)
+ IF (ips_i .LE. pmaxi_cp+1) THEN
+ DO ip=ips_i,pmaxloop_i
+ IF (ijs_i .LE. jmaxi_cp+1) THEN
+ DO ij=ijs_i,jmaxloop_i
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
+ DO iz = izs,ize
+ moments_i(ip,ij,ikx,iky,iz,:) = moments_i_cp(ip,ij,ikx,iky,iz)
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
+ ENDIF
+ ENDDO
+ ENDIF
! Deallocate checkpoint arrays
DEALLOCATE(moments_i_cp)
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 0d57eb336ae8322c3f325f0a5f265a715f5bee74..116f310eb3f5e23d2a90d18f702206b14a99d3b7 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='6d12f9c-dirty')
+parameter (VERSION='3b26011-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Tue Dec 14 10:47:37 CET 2021')
+parameter (EXECDATE='Fri Feb 25 14:24:41 CET 2022')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 0d57eb336ae8322c3f325f0a5f265a715f5bee74..116f310eb3f5e23d2a90d18f702206b14a99d3b7 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='6d12f9c-dirty')
+parameter (VERSION='3b26011-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Tue Dec 14 10:47:37 CET 2021')
+parameter (EXECDATE='Fri Feb 25 14:24:41 CET 2022')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index f56880aee24408625d08627b5cb56cd9a6081752..79f5dfbb6986931291ab374612848b0f5e5978e2 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -27,6 +27,7 @@ SUBROUTINE stepon
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
IF(KIN_E) CALL moments_eq_rhs_e
CALL moments_eq_rhs_i
+
! ---- step n -> n+1 transition
! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
CALL advance_time_level
@@ -37,12 +38,12 @@ SUBROUTINE stepon
CALL apply_closure_model
! Exchanges the ghosts values of N_n+1
CALL update_ghosts
- ! Update collision C_n+1 = C(N_n+1)
- CALL compute_TColl
! Update electrostatic potential phi_n = phi(N_n+1)
CALL poisson
! Update non adiabatic moments n -> n+1
CALL compute_nadiab_moments
+ ! Update collision C_n+1 = C(N_n+1)
+ CALL compute_TColl
! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
CALL compute_Sapj
! Store or cancel/maintain zonal modes artificially
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index d3e9d353f6a6a3698b0e1efce5c24683851163a4..568a56d7da05b2ca2539d4da693bd6b31dd2b27b 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -29,7 +29,7 @@ end
if 0
%% statistical transport averaging
-options.T = [1500 5000];
+options.T = [5000 6000];
fig = statistical_transport_averaging(data,options);
end
@@ -112,7 +112,7 @@ end
if 0
%% 1D spectrum
-options.TIME = 1000:3000;
+options.TIME = 5000:10:5200;
options.NORM = 1;
% options.NAME = '\phi';
% options.NAME = 'n_i';
@@ -147,7 +147,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 0.01:0.01:1.0;
-options.TIME =500:1:3000;
+options.TIME = 5000:1:5050;
options.NMA = 1;
options.NMODES = 20;
fig = mode_growth_meter(data,options);
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index aa87383b5a00f19bb9efbaa6606c92856d59b178..e3f8e09f2f843a2799989887c352fc9e29868b63 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -88,7 +88,7 @@ outfile ='';
%% nu = 0
-outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
@@ -138,8 +138,11 @@ outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e
% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
+% outfile = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1';
+
% outfile = 'predator_prey/SG_Kn_1.7_nu_0.01';
+outfile = 'ZF_damping_DK/SG_4x2_nu_0.1';
% else% Marconi results
% outfile ='';
% outfile ='';
diff --git a/wk/benchmark_HeLaZ_gene_transport_zpinch.m b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
index 0359ffb73741d5073bf2fc86533252e6d25714e3..59d756ebeb1a48be79b3f9fe1e8ab01957078656 100644
--- a/wk/benchmark_HeLaZ_gene_transport_zpinch.m
+++ b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
@@ -86,10 +86,10 @@ fig = figure; set(gcf,'Position',[250 250 600 300]);
% upward scan
KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.3e-1 1.3e+0 2.3e+0 3.6e+0 6.3e+0 9.4e+0 1.5e+1 2.3e+1 3.8e+1];
-LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.0e+0 3.5e+0 6.0e+0 4.8e+0 6.0e+0 2.9e+1];
+LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.3e+0 3.5e+0 6.0e+0 4.8e+0 6.0e+0 2.9e+1];
SGGK_200x32x05x03 = [0.0e-9 1.5e-2 1.0e-1 3.8e-1 8.4e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
-DGGK_200x32x05x03 = [2.0e-4 3.6e-3 1.3e-2 4.0e-2 1.2e-1 1.8e-1 2.0e-1 4.4e-1 6.6e-1 1.3e+1 3.2e+1];
-NOCO_200x32x05x03 = [0.0e+0 1.2e-2 2.6e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 0.0e+0 0.0e+0 0.0e+1 2.4e+1];
+DGGK_200x32x05x03 = [2.0e-4 3.6e-3 1.3e-2 4.0e-2 1.2e-1 1.8e-1 2.0e-1 2.2e+0 6.6e-1 1.3e+1 3.2e+1];
+NOCO_200x32x05x03 = [0.0e+0 1.2e-2 2.6e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 3.9e+0 6.9e+0 1.2e+1 2.4e+1];
% downward scan
% KN = [KN 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
% LDGK_200x32x05x03 = [LDGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
@@ -122,10 +122,11 @@ xlim([ 1.55 2.55]);
% saveas(fig,'/home/ahoffman/Dropbox/Applications/Overleaf/Paper 1/figures/coll_transport_benchmark.eps')
%% Burst study Kn = 1.7
-NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 1.0e-1];
-SGGK_200x32x05x03 = [0.0e+0 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 1.0e-1];
-LRGK_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e-2 0.0e-2 0.0e+0 0.0e+0];
+NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 6.0e-2 7.0e-2 8.0e-2 9.0e-2 1.0e-1];
+SGGK_200x32x05x03 = [1.0e-2 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 6.5e-2 7.5e-2 8.3e-2 8.8e-2 1.0e-1];
+SGGK_200x32x10x05 = [1.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+LRGK_200x32x05x03 = [0.0e+0 2.9e-2 0.0e+0 0.0e-2 0.0e-2 0.0e+0 0.0e+0];
% NOCO_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
-semilogy(NU,SGGK_200x32x05x03); hold on
+plot(NU,SGGK_200x32x05x03,'x'); hold on
grid on; xlabel('$\nu$'); ylabel('$\Gamma^\infty_x$');
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index 82ec97b90edcd0b88841ac038f2288e43e2cd8fb..e6d1265e21d30487fdd505030591353b4cd1de88 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -6,7 +6,7 @@ default_plots_options
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_N = 1.7; % Density gradient drive
K_T = 0.25*K_N; % Temperature '''
@@ -23,7 +23,7 @@ SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 300; % Maximal time unit
+TMAX = 100; % Maximal time unit
DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
@@ -37,8 +37,8 @@ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'SG';
-GKCO = 1; % gyrokinetic operator
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
@@ -72,11 +72,11 @@ CURVB = 1.0;
if 1
% Parameter scan over PJ
-PA = [4 6 8 10];
-JA = [2 3 4 5];
+% PA = [4 6 8 10];
+% JA = [2 3 4 5];
Nparam = numel(PA);
% Parameter scan over KN
-% PA = [4]; JA = [2];
+PA = [4]; JA = [2];
% PMAXE = PA(1); PMAXI = PA(1);
% JMAXE = JA(1); JMAXI = JA(1);
% KNA = 1.5:0.05:2.5;