diff --git a/Makefile b/Makefile
index 8d514dcc7d7318e0c35b61d39927c75427eb40af..c405b8e434025e986bb7984516a8f806f072c3d7 100644
--- a/Makefile
+++ b/Makefile
@@ -114,7 +114,7 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
- $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
$(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 65b13c91548e2b648e6f6c5281d70b7765f53ce8..32ac436755659dbdcee89b79ab50d3d9f5e3612d 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -7,4 +7,4 @@ Y_ = dataObjs(1).YData;
figure;
plot(X_,Y_);
-plot(X_(9000:12000)-X_(8000),Y_(9000:12000));
\ No newline at end of file
+plot(X_-X_(1),Y_);
\ No newline at end of file
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 1b0abf43828ebee86a1e80bd1bcf018a78873712..fe3e99259d224c5807cdfee2d41bc30c2e199d84 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -1,9 +1,12 @@
CO = h5readatt(filename,'/data/input','CO');
% K_N = h5readatt(filename,'/data/input','eta_n');
% K_T = h5readatt(filename,'/data/input','eta_T');
-K_N = h5readatt(filename,'/data/input','K_n');
-K_T = h5readatt(filename,'/data/input','K_T');
-K_E = h5readatt(filename,'/data/input','K_E');
+K_N = h5readatt(filename,'/data/input','K_n');
+K_T = h5readatt(filename,'/data/input','K_T');
+K_E = h5readatt(filename,'/data/input','K_E');
+Q0 = h5readatt(filename,'/data/input','q0');
+SHEAR = h5readatt(filename,'/data/input','shear');
+EPS = h5readatt(filename,'/data/input','eps');
PMAXI = h5readatt(filename,'/data/input','pmaxi');
JMAXI = h5readatt(filename,'/data/input','jmaxi');
PMAXE = h5readatt(filename,'/data/input','pmaxe');
diff --git a/matlab/load_results.m b/matlab/load_results.m
index 81bac625891aec776f3deb0eb56b96651937bc16..ad16253186baa5d97f83f08c35f7b7f72b5911c2 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -14,6 +14,8 @@ W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
+% KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
+KIN_E = 1;
if W_GAMMA
@@ -31,26 +33,36 @@ end
if W_NA00
[Ni00, Ts3D, dt3D] = load_3D_data(filename, 'Ni00');
+ if(KIN_E)
Ne00 = load_3D_data(filename, 'Ne00');
+ end
end
if W_NAPJ
[Nipj, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
+ if(KIN_E)
[Nepj ] = load_5D_data(filename, 'moments_e');
+ end
end
if W_SAPJ
[Sipj, Ts5D, dt5D] = load_5D_data(filename, 'Sipj');
+ if(KIN_E)
Sepj = load_5D_data(filename, 'Sepj');
+ end
end
if W_DENS
+ if(KIN_E)
[DENS_E, Ts3D, dt3D] = load_3D_data(filename, 'dens_e');
+ end
[DENS_I, Ts3D, dt3D] = load_3D_data(filename, 'dens_i');
end
if W_TEMP
+ if(KIN_E)
[TEMP_E, Ts3D, dt3D] = load_3D_data(filename, 'temp_e');
+ end
[TEMP_I, Ts3D, dt3D] = load_3D_data(filename, 'temp_i');
end
\ No newline at end of file
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index 434be93a4d18a9086d2fb1990652a50c414e89dc..1cccc49ff54b7d090c4c3ab127f40c78b5b99927 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -23,7 +23,6 @@ Lx = 2*pi/dkx; Ly = 2*pi/dky;
dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz;
x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
-z = dz * (1:Nz);
[Y_XY,X_XY] = meshgrid(y,x);
[Z_XZ,X_XZ] = meshgrid(z,x);
[Z_YZ,Y_YZ] = meshgrid(z,y);
@@ -153,45 +152,45 @@ for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
end
%% Compute primary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-[~,itmax] = min(abs(Ts3D-tmax));
-tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax);
-[~,its3D_lin] = min(abs(Ts3D-tstart));
-[~,ite3D_lin] = min(abs(Ts3D-tend));
-
-g_I = zeros(Nkx,Nky,Nz);
-for ikx = 1:Nkx
- for iky = 1:Nky
- for iz = 1:Nz
- [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
- end
- end
-end
-[gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
-kmax_I = abs(ky(ikmax_I));
-Bohm_transport = K_N*gmax_I/kmax_I^2;
+% disp('- growth rate')
+% % Find max value of transport (end of linear mode)
+% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
+% [~,itmax] = min(abs(Ts3D-tmax));
+% tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax);
+% [~,its3D_lin] = min(abs(Ts3D-tstart));
+% [~,ite3D_lin] = min(abs(Ts3D-tend));
+%
+% g_I = zeros(Nkx,Nky,Nz);
+% for ikx = 1:Nkx
+% for iky = 1:Nky
+% for iz = 1:Nz
+% [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+% end
+% end
+% end
+% [gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
+% kmax_I = abs(ky(ikmax_I));
+% Bohm_transport = K_N*gmax_I/kmax_I^2;
%% Compute secondary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-% [~,itmax] = min(abs(Ts2D-tmax));
-% tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
-[~,its3D_lin] = min(abs(Ts3D-tstart));
-[~,ite3D_lin] = min(abs(Ts3D-tend));
-
-g_II = zeros(Nkx,Nky);
-for ikx = 1:Nkx
- for iky = 1
- for iz = 1:Nz
- [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
- end
- end
-end
-[gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
-kmax_II = abs(kx(ikmax_II));
+% disp('- growth rate')
+% % Find max value of transport (end of linear mode)
+% % [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
+% % [~,itmax] = min(abs(Ts2D-tmax));
+% % tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
+% [~,its3D_lin] = min(abs(Ts3D-tstart));
+% [~,ite3D_lin] = min(abs(Ts3D-tend));
+%
+% g_II = zeros(Nkx,Nky);
+% for ikx = 1:Nkx
+% for iky = 1
+% for iz = 1:Nz
+% [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
+% end
+% end
+% end
+% [gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
+% kmax_II = abs(kx(ikmax_II));
%% zonal vs nonzonal energies for phi(t)
Ephi_Z = zeros(1,Ns3D);
@@ -210,7 +209,7 @@ for it = 1:numel(Ts3D)
% Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2;
Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2))));
% Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it);
- high_k_phi(it) = squeeze(abs(PHI(18,18,1,it)).^2); % kperp sqrt(2)
+% high_k_phi(it) = squeeze(abs(PHI(18,18,1,it)).^2); % kperp sqrt(2)
% high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5
end
diff --git a/matlab/setup.m b/matlab/setup.m
index e272a4f33cd119fc90ebf712e7033925291f8434..a9a1f15ba5d1f0e8f5468264b2ea9bf42c07db47 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -5,20 +5,22 @@ GRID.pmaxe = PMAXE; % Electron Hermite moments
GRID.jmaxe = JMAXE; % Electron Laguerre moments
GRID.pmaxi = PMAXI; % Ion Hermite moments
GRID.jmaxi = JMAXI; % Ion Laguerre moments
-GRID.Nx = Nx; % x grid resolution
-GRID.Lx = Lx; % x length
-GRID.Ny = Ny; % y ''
-GRID.Ly = Ly; % y ''
-GRID.Nz = Nz; % z resolution
-GRID.q0 = q0; % q factor
-GRID.shear = shear; % shear
-GRID.eps = eps; % inverse aspect ratio
+GRID.Nx = NX; % x grid resolution
+GRID.Lx = LX; % x length
+GRID.Ny = NY; % y ''
+GRID.Ly = LY; % y ''
+GRID.Nz = NZ; % z resolution
+GRID.q0 = Q0; % q factor
+GRID.shear = SHEAR; % shear
+GRID.eps = EPS; % inverse aspect ratio
% Model parameters
MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
+MODEL.KIN_E = KIN_E;
+if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
MODEL.mu = MU;
MODEL.mu_p = MU_P;
MODEL.mu_j = MU_J;
@@ -98,28 +100,30 @@ coll_ = sprintf(coll_,NU);
% nonlinear
lin_ = [];
if ~NON_LIN; lin_ = '_lin'; end
+adiabe_ = [];
+if ~KIN_E; adiabe_ = '_adiabe'; end
% resolution and boxsize
res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)];
-if (Lx ~= Ly)
- geo_ = ['_Lx_',num2str(Lx),'_Ly_',num2str(Ly)];
+if (LX ~= LY)
+ geo_ = ['_Lx_',num2str(LX),'_Ly_',num2str(LY)];
else
- geo_ = ['_L_',num2str(Lx)];
+ geo_ = ['_L_',num2str(LX)];
end
-if (Nz > 1) %3D case
- res_ = [res_,'x',num2str(GRID.Nz)];
- if abs(q0) > 0
- geo_ = [geo_,'_q0_',num2str(q0)];
+if (NZ > 1) %3D case
+ res_ = [res_,'x',num2str(NZ)];
+ if abs(Q0) > 0
+ geo_ = [geo_,'_q0_',num2str(Q0)];
end
- if abs(eps) > 0
- geo_ = [geo_,'_e_',num2str(eps)];
+ if abs(EPS) > 0
+ geo_ = [geo_,'_e_',num2str(EPS)];
end
- if abs(shear) > 0
- geo_ = [geo_,'_s_',num2str(shear)];
+ if abs(SHEAR) > 0
+ geo_ = [geo_,'_s_',num2str(SHEAR)];
end
end
% put everything together in the param character chain
u_ = '_'; % underscore variable
-PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_];
+PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_];
BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/'];
BASIC.PARAMS = PARAMS;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 57873487a7619ca979d19fb1059acbbb446f92cb..84259b0f77cca3641911abaf13592dc1cf25f820 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -49,6 +49,7 @@ fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']);
fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']);
fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']);
fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']);
+fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']);
fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 4f0aa5b9a94e7e8dfcdfed4c0df3e4a14871c5ee..5e8c87d622973eb76990cfba6658d4a3711d49c4 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -20,20 +20,22 @@ CONTAINS
USE time_integration
USE grid
use prec_const
- USE model, ONLY: CLOS
+ USE model, ONLY: CLOS, KIN_E
use fields, ONLY: moments_e, moments_i
use array, ONLY: moments_rhs_e, moments_rhs_i
IMPLICIT NONE
INTEGER :: p_int, j_int
CALL cpu_time(t0_adv_field)
- DO ip=ips_e,ipe_e
- p_int = parray_e(ip)
- DO ij=ijs_e,ije_e
- IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
- CALL advance_field(moments_e(ip,ij,:,:,:,:), moments_rhs_e(ip,ij,:,:,:,:))
+ IF(KIN_E) THEN
+ DO ip=ips_e,ipe_e
+ p_int = parray_e(ip)
+ DO ij=ijs_e,ije_e
+ IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
+ CALL advance_field(moments_e(ip,ij,:,:,:,:), moments_rhs_e(ip,ij,:,:,:,:))
+ ENDDO
ENDDO
- ENDDO
+ ENDIF
DO ip=ips_i,ipe_i
p_int = parray_i(ip)
DO ij=ijs_i,ije_i
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 9c14fcd3f0c40037c45549a3465b91b046eef58c..668a95bcca7c4fa0b3467c4496db5de57af75eb2 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -40,8 +40,6 @@ MODULE array
! Geoemtrical operators
! Curvature
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z
- ! kperp array
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, PUBLIC :: kparray ! dimensions: kx, ky, z
! Jacobian
REAL(dp), DIMENSION(:), ALLOCATABLE :: Jacobian ! dimensions: z
! Metric
@@ -62,6 +60,8 @@ MODULE array
! Kernel function evaluation (ij,ikx,iky,iz)
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_e
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_i
+ ! Poisson operator (ikx,iky,iz)
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
!! Diagnostics
! Gyrocenter density for electron and ions (ikx,iky,iz)
diff --git a/src/auxval.F90 b/src/auxval.F90
index ea051ab30482a2960c1baeb8a95b4d06e63a8b9a..8e642f6f8b7e744540429c938332ace97c5aedd3 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -21,7 +21,7 @@ subroutine auxval
ENDIF
! Init the grids
CALL set_pgrid ! parallel kin (MPI distributed)
-
+
CALL set_jgrid ! perp kin
CALL set_kxgrid ! radial modes (MPI distributed by FFTW)
@@ -38,6 +38,8 @@ subroutine auxval
CALL evaluate_kernels ! precompute the kernels
+ CALL evaluate_poisson_op ! precompute the kernels
+
IF ( NON_LIN ) THEN;
CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs
ENDIF
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index d3c48200e77ac6c39be794af4a903c5127e1c698..4b124fc32c7c2abfe9e397e69d27566bd26d9012 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -1,7 +1,7 @@
module closure
! Contains the routines to define closures
USE basic
-USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, nu
+USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, nu, KIN_E
USE grid
USE array, ONLY: kernel_e, kernel_i
USE fields, ONLY: moments_e, moments_i
@@ -31,12 +31,14 @@ SUBROUTINE apply_closure_model
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
+ IF(KIN_E) THEN
DO ip = ipsg_e,ipeg_e
DO ij = ijsg_e,ijeg_e
IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) &
moments_e(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
ENDDO
+ ENDIF
DO ip = ipsg_i,ipeg_i
DO ij = ijsg_i,ijeg_i
IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) &
@@ -65,6 +67,7 @@ SUBROUTINE ghosts_truncation
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
+ IF(KIN_E) THEN
DO ip = ipsg_e,ipeg_e
moments_e(ip,ijsg_e,ikx,iky,iz,updatetlevel) = 0._dp
moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp
@@ -79,7 +82,7 @@ SUBROUTINE ghosts_truncation
ENDDO
kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
-
+ ENDIF
DO ip = ipsg_i,ipeg_i
moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 4a0a291113656f220d4d5be7e36e5549222a8b12..3d2bd674d1358deaf2d320569451d7bfa0430364 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -67,8 +67,8 @@ CONTAINS
!******************************************************************************!
SUBROUTINE DoughertyGK_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
USE fields, ONLY: moments_e, phi
- USE grid, ONLY: parray_e, jarray_e, kxarray, kyarray, Jmaxe, ip0_e, ip1_e, ip2_e
- USE array, ONLY: kernel_e, kparray
+ 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
@@ -173,8 +173,8 @@ CONTAINS
!******************************************************************************!
SUBROUTINE DoughertyGK_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
USE fields, ONLY: moments_i, phi
- USE grid, ONLY: parray_i, jarray_i, kxarray, kyarray, Jmaxi, ip0_i, ip1_i, ip2_i
- USE array, ONLY: kernel_i, kparray
+ USE grid, ONLY: parray_i, jarray_i, kxarray, kyarray, kparray, Jmaxi, ip0_i, ip1_i, ip2_i
+ USE array, ONLY: kernel_i
USE basic
USE model, ONLY: sigmai2_taui_o2, qi_taui, nu_i, sqrt_sigmai2_taui_o2
USE time_integration, ONLY : updatetlevel
@@ -303,6 +303,7 @@ CONTAINS
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
+ IF(KIN_E) THEN
! Electrons
DO ij = 1,Jmaxe+1
! Loop over all p to compute sub collision term
@@ -325,6 +326,7 @@ CONTAINS
TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip)
ENDDO
ENDDO
+ ENDIF
! Ions
DO ij = 1,Jmaxi+1
DO ip = 1,total_np_i
@@ -419,7 +421,7 @@ CONTAINS
USE basic
USE time_integration, ONLY : updatetlevel
USE utility
- USE model, ONLY: CO, nu_i, nu_ie, CLOS
+ USE model, ONLY: CO, nu_i, nu_ie, CLOS, KIN_E
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -440,12 +442,16 @@ CONTAINS
jloopii: DO ij2 = ijs_i,ije_i
j2_int = jarray_i(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
+ ! Ion-ion collision
TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
- *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C) &
- +nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C))
+ * nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C)
+ IF(KIN_E) & ! Ion-electron collision test
+ TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ * nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C)
ENDDO jloopii
ENDDO ploopii
+ IF(KIN_E) THEN ! Ion-electron collision field
ploopie: DO ip2 = ips_e,ipe_e ! sum the ion-electron field terms
p2_int = parray_e(ip2)
jloopie: DO ij2 = ijs_e,ije_e
@@ -455,6 +461,7 @@ CONTAINS
*(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikx_C, iky_C))
ENDDO jloopie
ENDDO ploopie
+ ENDIF
END SUBROUTINE apply_COSOlver_mat_i
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 967e1039db9e1f93187521bd8857382a73b38297..bfc834960b7c5604d0729a9a6df3e0c3e20ca0e7 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -25,8 +25,9 @@ SUBROUTINE compute_Sapj
! Execution time start
CALL cpu_time(t0_Sapj)
-zloop: DO iz = izs,ize
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
+ IF(KIN_E) THEN
+ zloope: DO iz = izs,ize
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
p_int = parray_e(ip)
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
@@ -118,9 +119,12 @@ zloop: DO iz = izs,ize
ENDIF GF_CLOSURE_e
ENDDO jloope
ENDDO ploope
+ENDDO zloope
+ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
+zloopi: DO iz = izs,ize
ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
! we check if poly degree is even (eq to index is odd) to spare computation
@@ -214,11 +218,10 @@ zloop: DO iz = izs,ize
ENDIF GF_CLOSURE_i
ENDDO jloopi
ENDDO ploopi
+ENDDO zloopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Execution time END
CALL cpu_time(t1_Sapj)
tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
-
-ENDDO zloop
END SUBROUTINE compute_Sapj
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 4ab9b96142955695cff335cbf24cdcff602cc44a..1ba38d81a0e333fa1b5757ea43fc87411a70157b 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -109,16 +109,19 @@ SUBROUTINE diagnose(kstep)
IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
IF (write_Na00) THEN
+ IF(KIN_E)&
CALL creatg(fidres, "/data/var3d/Ne00", "Ne00")
CALL creatg(fidres, "/data/var3d/Ni00", "Ni00")
ENDIF
IF (write_dens) THEN
+ IF(KIN_E)&
CALL creatg(fidres, "/data/var3d/dens_e", "dens_e")
CALL creatg(fidres, "/data/var3d/dens_i", "dens_i")
ENDIF
IF (write_temp) THEN
+ IF(KIN_E)&
CALL creatg(fidres, "/data/var3d/temp_e", "temp_e")
CALL creatg(fidres, "/data/var3d/temp_i", "temp_i")
ENDIF
@@ -136,11 +139,13 @@ SUBROUTINE diagnose(kstep)
CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
IF (write_Napj) THEN
+ IF(KIN_E)&
CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
ENDIF
IF (write_Sapj) THEN
+ IF(KIN_E)&
CALL creatg(fidres, "/data/var5d/moments_e", "Sipj")
CALL creatg(fidres, "/data/var5d/moments_i", "Sepj")
ENDIF
@@ -291,6 +296,7 @@ SUBROUTINE diagnose_0d
USE diagnostics_par
USE prec_const
USE processing
+ USE model, ONLY: KIN_E
IMPLICIT NONE
! Time measurement data
@@ -335,6 +341,7 @@ SUBROUTINE diagnose_2d
USE diagnostics_par
USE prec_const
USE processing
+ USE model, ONLY: KIN_E
IMPLICIT NONE
@@ -389,6 +396,7 @@ SUBROUTINE diagnose_3d
USE diagnostics_par
USE prec_const
USE processing
+ USE model, ONLY: KIN_E
IMPLICIT NONE
@@ -403,26 +411,30 @@ SUBROUTINE diagnose_3d
IF (write_phi) CALL write_field3d(phi (:,:,:), 'phi')
IF (write_Na00) THEN
- IF ( (ips_e .EQ. 1) .AND. (ips_i .EQ. 1) ) THEN
+ IF(KIN_E)THEN
+ IF (ips_e .EQ. 1) THEN
Ne00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_e(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
- Ni00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_i(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
ENDIF
-
CALL manual_3D_bcast(Ne00(ikxs:ikxe,ikys:ikye,izs:ize))
CALL write_field3d(Ne00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ne00')
-
+ ENDIF
+ IF (ips_i .EQ. 1) THEN
+ Ni00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_i(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ ENDIF
CALL manual_3D_bcast(Ni00(ikxs:ikxe,ikys:ikye,izs:ize))
CALL write_field3d(Ni00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ni00')
ENDIF
IF (write_dens) THEN
CALL compute_density
+ IF(KIN_E)&
CALL write_field3d(dens_e(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_e')
CALL write_field3d(dens_i(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_i')
ENDIF
IF (write_temp) THEN
CALL compute_temperature
+ IF(KIN_E)&
CALL write_field3d(temp_e(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_e')
CALL write_field3d(temp_i(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_i')
ENDIF
@@ -463,6 +475,7 @@ SUBROUTINE diagnose_5d
USE time_integration
USE diagnostics_par
USE prec_const
+ USE model, ONLY: KIN_E
IMPLICIT NONE
CALL append(fidres, "/data/var5d/time", time,ionode=0)
@@ -472,11 +485,13 @@ SUBROUTINE diagnose_5d
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
IF (write_Napj) THEN
+ IF(KIN_E)&
CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,:,updatetlevel), 'moments_e')
CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,:,updatetlevel), 'moments_i')
ENDIF
IF (write_Sapj) THEN
+ IF(KIN_E)&
CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:,:), 'Sepj')
CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:,:), 'Sipj')
ENDIF
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index be73c0a59664258d56062cb0ffb0a20b566b18be..2ef6d63d89deb20184000084cd792aa9a8d09b75 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -8,9 +8,11 @@ module geometry
use array
use fields
use basic
+ use utility, ONLY: simpson_rule_z
implicit none
public
+ COMPLEX(dp), PROTECTED :: iInt_Jacobian
contains
@@ -18,6 +20,7 @@ contains
! evalute metrix, elements, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
+ COMPLEX(dp), DIMENSION(izs:ize) :: integrant
!
IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
@@ -41,6 +44,11 @@ contains
ENDDO
ENDDO
ENDDO
+ !
+ ! Compute the inverse z integrated Jacobian (useful for flux averaging)
+ integrant = Jacobian ! Convert into complex array
+ CALL simpson_rule_z(integrant,iInt_Jacobian)
+ iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index c3a4e82d92f404b8b1980a4e18ab23b1a9f56c50..fa0757f9b15a1a7b26d92027bbd3934aeee189e8 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -3,6 +3,7 @@ USE basic
USE fields, ONLY : moments_e, moments_i
USE grid
USE time_integration
+USE model, ONLY : KIN_E
IMPLICIT NONE
@@ -17,7 +18,7 @@ SUBROUTINE update_ghosts
IF (num_procs_p .GT. 1) THEN ! Do it only if we share the p
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
- CALL update_ghosts_p_e
+ IF(KIN_E) CALL update_ghosts_p_e
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
CALL update_ghosts_p_i
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 9a6d695b3bea3117884ab29b2f7fe56fa4de04f1..a2736db8c102e24b5204804c3cd7216f61cc8157 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -57,6 +57,7 @@ MODULE grid
! Grids containing position in fourier space
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kxarray, kxarray_full
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kyarray, kyarray_full
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, PUBLIC :: kparray
REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max!, kp_max
REAL(dp), PUBLIC, PROTECTED :: local_kxmax, local_kymax
INTEGER, PUBLIC, PROTECTED :: ikxs, ikxe, ikys, ikye!, ikps, ikpe
diff --git a/src/inital.F90 b/src/inital.F90
index 3bac2e8902e2a8751451f18e28b653d796089603..0627b10e52124356cc6bf32e09299a9548ddd6b7 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -274,7 +274,7 @@ SUBROUTINE init_phi
IMPLICIT NONE
REAL(dp) :: noise
- REAL(dp) :: kx, ky, sigma, gain, ky_shift
+ REAL(dp) :: kx, ky, kp, sigma, gain, ky_shift
INTEGER, DIMENSION(12) :: iseedarr
! Seed random number generator
@@ -286,8 +286,9 @@ SUBROUTINE init_phi
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
+ kp = kparray(ikx,iky,iz)
CALL RANDOM_NUMBER(noise)
- phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
+ phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))*EXP(-0.1*kp**2)!*AA_x(ikx)*AA_y(iky)
ENDDO
END DO
END DO
diff --git a/src/memory.F90 b/src/memory.F90
index e6dc6a2abb935e9721364d5c2e9247b384d6722e..0cfc3203f5d87ed3d4fc474e0ca2afd8a09509d9 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -6,70 +6,52 @@ SUBROUTINE memory
USE fields
USE grid
USE time_integration
- USE model, ONLY: CO, NON_LIN
+ USE model, ONLY: CO, NON_LIN, KIN_E
USE prec_const
IMPLICIT NONE
- !___________________ 2D ARRAYS __________________________
+
! Electrostatic potential
CALL allocate_array(phi, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izs,ize)
+
+ !Electrons arrays
+ IF(KIN_E) THEN
+ CALL allocate_array( Ne00, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( dens_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( xnepj, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( xnepp2j, ips_e,ipe_e)
+ CALL allocate_array( xnepp1j, ips_e,ipe_e)
+ CALL allocate_array( xnepm1j, ips_e,ipe_e)
+ CALL allocate_array( xnepm2j, ips_e,ipe_e)
+ CALL allocate_array( xnepjp1, ijs_e,ije_e)
+ CALL allocate_array( xnepjm1, ijs_e,ije_e)
+ CALL allocate_array( ynepp1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepm1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepp1jm1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1jp1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
+ ENDIF
- !! Diagnostics arrays
- ! Gyrocenter density *for 2D output*
- CALL allocate_array(Ne00, ikxs,ikxe, ikys,ikye, izs,ize)
+ !Ions arrays
CALL allocate_array(Ni00, ikxs,ikxe, ikys,ikye, izs,ize)
- ! particle density *for 2D output*
- CALL allocate_array(dens_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(dens_i, ikxs,ikxe, ikys,ikye, izs,ize)
- ! particle temperature *for 2D output*
- CALL allocate_array(temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
-
- !___________________ 4D ARRAYS __________________________
- !! FLR kernels functions
- ! Kernel evaluation from j= -1 to jmax+1 for truncation
- CALL allocate_array(Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
-
- !___________________ 5D ARRAYS __________________________
- ! Moments with ghost degrees for p+2 p-2, j+1, j-1 truncations
- CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
-
- ! Moments right-hand-side (contains linear part of hierarchy)
- CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
-
- ! Non linear terms and dnjs table
- CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
-
- ! Collision term
- CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i , ikxs,ikxe, ikys,ikye, izs,ize)
-
- ! Non linear terms and dnjs table
- CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
-
- ! Linear coeff for moments rhs
- ! electrons
- CALL allocate_array( xnepj, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( xnepp2j, ips_e,ipe_e)
- CALL allocate_array( xnepp1j, ips_e,ipe_e)
- CALL allocate_array( xnepm1j, ips_e,ipe_e)
- CALL allocate_array( xnepm2j, ips_e,ipe_e)
- CALL allocate_array( xnepjp1, ijs_e,ije_e)
- CALL allocate_array( xnepjm1, ijs_e,ije_e)
- CALL allocate_array( ynepp1j, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( ynepm1j, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( ynepp1jm1, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( ynepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( zNepm1j, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( zNepm1jp1, ips_e,ipe_e, ijs_e,ije_e)
- CALL allocate_array( zNepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
- ! ions
+ CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( xnipj, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xnipp2j, ips_i,ipe_i)
CALL allocate_array( xnipp1j, ips_i,ipe_i)
@@ -84,7 +66,11 @@ SUBROUTINE memory
CALL allocate_array( zNipm1j, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( zNipm1jp1, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( zNipm1jm1, ips_i,ipe_i, ijs_i,ije_i)
- ! elect. pot.
+
+ ! Non linear terms and dnjs table
+ CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
+
+ ! elect. pot. linear terms
CALL allocate_array( xphij, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijp1, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
@@ -110,19 +96,23 @@ SUBROUTINE memory
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
IF (CO .LT. -1) THEN !DK collision matrix (same for every k)
+ IF (KIN_E) THEN
CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
- CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ ENDIF
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
ELSEIF (CO .GT. 1) THEN !GK collision matrices (one for each kperp)
+ IF (KIN_E) THEN
CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye)
- CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye)
CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye)
CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
+ ENDIF
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye)
ENDIF
END SUBROUTINE memory
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 10bde0bbd344f6c0aec6e96aad3f33127118385b..dac8b664894f6833d129a574e4401b9f91389abd 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -10,6 +10,7 @@ MODULE model
INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method
INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model
LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term
+ LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise)
REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! spatial Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
@@ -38,7 +39,7 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: sigmae2_taue_o2 ! factor of the Kernel argument
REAL(dp), PUBLIC, PROTECTED :: sigmai2_taui_o2 !
REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmae2_taue_o2 ! factor of the Kernel argument
- REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmai2_taui_o2
+ REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmai2_taui_o2
REAL(dp), PUBLIC, PROTECTED :: nu_e, nu_i ! electron-ion, ion-ion collision frequency
REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie ! e-e, i-e coll. frequ.
REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum
@@ -53,7 +54,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CO, CLOS, NL_CLOS, KERN, NON_LIN, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CO, CLOS, NL_CLOS, KERN, NON_LIN, KIN_E, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, K_n, K_T, K_E, GradB, CurvB, lambdaD
READ(lu_in,model_par)
@@ -102,6 +103,7 @@ CONTAINS
CALL attach(fidres, TRIM(str), "CLOS", CLOS)
CALL attach(fidres, TRIM(str), "KERN", KERN)
CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN)
+ CALL attach(fidres, TRIM(str), "KIN_E", KIN_E)
CALL attach(fidres, TRIM(str), "nu", nu)
CALL attach(fidres, TRIM(str), "mu", mu)
CALL attach(fidres, TRIM(str), "tau_e", tau_e)
@@ -110,9 +112,9 @@ CONTAINS
CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
CALL attach(fidres, TRIM(str), "q_e", q_e)
CALL attach(fidres, TRIM(str), "q_i", q_i)
- CALL attach(fidres, TRIM(str), "K_n", K_n)
- CALL attach(fidres, TRIM(str), "K_T", K_T)
- CALL attach(fidres, TRIM(str), "K_E", K_E)
+ CALL attach(fidres, TRIM(str), "K_n", K_n)
+ CALL attach(fidres, TRIM(str), "K_T", K_T)
+ CALL attach(fidres, TRIM(str), "K_E", K_E)
CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
END SUBROUTINE model_outputinputs
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 1660a18484201e08ca799f119ba0456d3f7b17b8..34a424de5adbfc3776ac3d27693d5de0888c58c1 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -68,6 +68,8 @@ SUBROUTINE moments_eq_rhs_e
! term propto n_e^{p,j-1}
Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,ikx,iky,iz)
! Parallel dynamic
+ Tpare = 0._dp; Tmir = 0._dp
+ IF(Nz .GT. 1) THEN
! ddz derivative for Landau damping term
Tpare = xnepp1j(ip) * &
( onetwelfth*nadiab_moments_e(ip+1,ij,ikx,iky,izm2)&
@@ -90,7 +92,7 @@ SUBROUTINE moments_eq_rhs_e
UNepm1jm1 = zNepm1jm1(ip,ij) * nadiab_moments_e(ip-1,ij-1,ikx,iky,iz)
Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1
-
+ ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij,ikx,iky,iz) &
@@ -213,6 +215,8 @@ SUBROUTINE moments_eq_rhs_i
! term propto n_e^{p,j-1}
Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip,ij-1,ikx,iky,iz)
! Parallel dynamic
+ Tpari = 0._dp; Tmir = 0._dp
+ IF(Nz .GT. 1) THEN
! term propto N_i^{p,j+1}, centered FDF
Tpari = xnipp1j(ip) * &
( onetwelfth*nadiab_moments_i(ip+1,ij,ikx,iky,izm2)&
@@ -235,7 +239,7 @@ SUBROUTINE moments_eq_rhs_i
Unipm1jm1 = znipm1jm1(ip,ij) * nadiab_moments_i(ip-1,ij-1,ikx,iky,iz)
Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1
-
+ ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij,ikx,iky,iz) &
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index edd297b22578b08c921d1cbd39c9f77efcd1b944..06235c2fae37f67ebcef16cc6af8f5708e6655f0 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -6,16 +6,11 @@ MODULE numerics
USE coeff
implicit none
- PUBLIC :: compute_derivatives, build_dnjs_table, evaluate_kernels, compute_lin_coeff
+ PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, compute_lin_coeff
PUBLIC :: wipe_turbulence, wipe_zonalflow
CONTAINS
-! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
-SUBROUTINE compute_derivatives
-
-END SUBROUTINE compute_derivatives
-
!******************************************************************************!
!!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
!******************************************************************************!
@@ -54,9 +49,10 @@ END SUBROUTINE build_dnjs_table
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
- USE array, Only : kernel_e, kernel_i, kparray
+ USE array, Only : kernel_e, kernel_i
USE grid
- USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
+ USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, &
+ lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2, KIN_E
IMPLICIT NONE
INTEGER :: j_int
REAL(dp) :: j_dp, y_, kp2_, kx_, ky_
@@ -65,12 +61,14 @@ DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
!!!!! Electron kernels !!!!!
+ IF(KIN_E) THEN
DO ij = ijsg_e, ijeg_e
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp)
y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz)**2
kernel_e(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
+ ENDIF
!!!!! Ion kernels !!!!!
DO ij = ijsg_i, ijeg_i
j_int = jarray_i(ij)
@@ -85,10 +83,55 @@ ENDDO
END SUBROUTINE evaluate_kernels
!******************************************************************************!
+!******************************************************************************!
+!!!!!!! Evaluate polarisation operator for Poisson equation
+!******************************************************************************!
+SUBROUTINE evaluate_poisson_op
+ USE basic
+ USE array, Only : kernel_e, kernel_i, inv_poisson_op
+ USE grid
+ USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E
+ IMPLICIT NONE
+ REAL(dp) :: pol_i, pol_e ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
+ INTEGER :: ini,ine
+
+ kxloop: DO ikx = ikxs,ikxe
+ kyloop: DO iky = ikys,ikye
+ zloop: DO iz = izs,ize
+ IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
+ inv_poisson_op(ikx, iky, iz) = 0._dp
+ ELSE
+ !!!!!!!!!!!!!!!!! Ion contribution
+ ! loop over n only if the max polynomial degree
+ pol_i = 0._dp
+ DO ini=1,jmaxi+1
+ pol_i = pol_i + qi2_taui*kernel_i(ini,ikx,iky,iz)**2 ! ... sum recursively ...
+ END DO
+ !!!!!!!!!!!!! Electron contribution\
+ pol_e = 0._dp
+ ! Kinetic model
+ IF (KIN_E) THEN
+ ! loop over n only if the max polynomial degree
+ DO ine=1,jmaxe+1 ! ine = n+1
+ pol_e = pol_e + qe2_taue*kernel_e(ine,ikx,iky,iz)**2 ! ... sum recursively ...
+ END DO
+ ! Adiabatic model
+ ELSE
+ pol_e = 1._dp - qe2_taue
+ ENDIF
+ inv_poisson_op(ikx, iky, iz) = 1._dp/(qe2_taue + qi2_taui - pol_i - pol_e)
+ ENDIF
+ END DO zloop
+ END DO kyloop
+ END DO kxloop
+
+END SUBROUTINE evaluate_poisson_op
+!******************************************************************************!
+
SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
- K_T, K_n, CurvB, GradB
+ K_T, K_n, CurvB, GradB, KIN_E
USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
@@ -97,6 +140,7 @@ SUBROUTINE compute_lin_coeff
REAL(dp) :: p_dp, j_dp
REAL(dp) :: kx, ky, z
!! Electrons linear coefficients for moment RHS !!!!!!!!!!
+ IF(KIN_E)THEN
DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
@@ -133,6 +177,7 @@ SUBROUTINE compute_lin_coeff
xnepjp1(ij) = -taue_qe * GradB * (j_dp + 1._dp)
xnepjm1(ij) = -taue_qe * GradB * j_dp
ENDDO
+ ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ions linear coefficients for moment RHS !!!!!!!!!!
DO ip = ips_i, ipe_i
diff --git a/src/poisson.F90 b/src/poisson.F90
index 486dc0cfda9f11ec3510a115615c0400a5da289d..8cc242ca37e3e6b4d6f0bc2bf650436554b0bb81 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -7,57 +7,61 @@ SUBROUTINE poisson
USE fields
USE grid
USE utility
- use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD
-
+ use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD, KIN_E
+ USE processing, ONLY : compute_density
USE prec_const
+ USE geometry, ONLY : iInt_Jacobian, Jacobian
IMPLICIT NONE
INTEGER :: ini,ine, i_, root_bcast
- REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build
- REAL(dp) :: polarisation ! sum_a(Z_a^2/tau_a (1-sum_n kernel_na^2))
- COMPLEX(dp) :: q_density ! charge density sum_a q_a n_a
- REAL(dp) :: gammaD
- COMPLEX(dp) :: gammaD_phi
+ COMPLEX(dp) :: fa_phi, intf_ ! current flux averaged phi
INTEGER :: count !! mpi integer to broadcast the electric potential at the end
COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye)
+ COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e ! charge density q_a n_a
+ COMPLEX(dp), DIMENSION(izs:ize) :: integrant ! charge density q_a n_a
!! Poisson can be solved only for process containing ip=1
IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
-
! Execution time start
CALL cpu_time(t0_poisson)
kxloop: DO ikx = ikxs,ikxe
kyloop: DO iky = ikys,ikye
- zloop: DO iz = izs,ize
+ !!!! Compute ion gyro charge density
+ rho_i = 0._dp
+ DO ini=1,jmaxi+1
+ rho_i = rho_i &
+ +q_i*kernel_i(ini,ikx,iky,:)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)
+ END DO
- q_density = 0._dp
- polarisation = 0._dp
- !!!!!!!!!!!!! Electron contribution
- ! loop over n only if the max polynomial degree
- DO ine=1,jmaxe+1 ! ine = n+1
- Kne = kernel_e(ine,ikx,iky,iz)
- q_density = q_density + q_e*Kne*moments_e(ip0_e, ine, ikx, iky, iz, updatetlevel)
- polarisation = polarisation + qe2_taue*Kne**2 ! ... sum recursively ...
+ !!!! Compute electron gyro charge density
+ rho_e = 0._dp
+ IF (KIN_E) THEN ! Kinetic electrons
+ DO ine=1,jmaxe+1
+ rho_e = rho_e &
+ +q_e*kernel_e(ine,ikx,iky,:)*moments_e(ip0_e,ine, ikx,iky,:,updatetlevel)
END DO
+ ELSE ! Adiabatic electrons
+ ! Adiabatic charge density (linked to flux averaged phi)
+ fa_phi = 0._dp
+ IF(kyarray(iky).EQ.0._dp) THEN
+ DO ini=1,jmaxi+1
+ rho_e(:) = Jacobian(:)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)&
+ *kernel_i(ini,ikx,iky,:)*(inv_poisson_op(ikx,iky,:)-1._dp)
+ call simpson_rule_z(rho_e,intf_)
+ fa_phi = fa_phi + intf_
+ ENDDO
+ rho_e = fa_phi*iInt_Jacobian !Normalize by 1/int(Jxyz)dz
+ ENDIF
+ ENDIF
- !!!!!!!!!!!!!!!!! Ions contribution
- ! loop over n only if the max polynomial degree
- DO ini=1,jmaxi+1
- Kni = kernel_i(ini,ikx,iky,iz)
- q_density = q_density + q_i*Kni*moments_i(ip0_i, ini, ikx, iky, iz, updatetlevel)
- polarisation = polarisation + qi2_taui*Kni**2 ! ... sum recursively ...
- END DO
-
- !!!!!!!!!!!!!!! Assembling the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
- phi(ikx, iky, iz) = q_density/(qe2_taue + qi2_taui - polarisation)
-
- END DO zloop
+ !!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
+ phi(ikx, iky, :) = (rho_e + rho_i)*inv_poisson_op(ikx,iky,:)
END DO kyloop
END DO kxloop
! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) phi(ikx_0,iky_0,izs:ize) = 0._dp
+ IF (contains_kx0 .AND. contains_ky0) phi(ikx_0,iky_0,:) = 0._dp
ENDIF
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 768cc94315754a43bec2948f330672168099dcfe..2866a44d0fad2ba3d056b1c836f2e59fa223a0b6 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -71,7 +71,7 @@ SUBROUTINE compute_radial_heatflux
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i
+ USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
IMPLICIT NONE
COMPLEX(dp) :: hflux_local
REAL(dp) :: ky_, buffer(1:2), j_dp
@@ -91,6 +91,7 @@ SUBROUTINE compute_radial_heatflux
* (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+ SQRT2/3._dp * kernel_i(ij,ikx,iky,iz) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel))
ENDDO
+ IF(KIN_E) THEN
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
hflux_local = hflux_local - imagu*ky_*CONJG(phi(ikx,iky,iz))&
@@ -101,6 +102,7 @@ SUBROUTINE compute_radial_heatflux
+q_e/tau_e * kernel_e( ij,ikx,iky,iz) * phi(ikx,iky,iz)) &
+SQRT2/3._dp * kernel_e(ij,ikx,iky,iz) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel))
ENDDO
+ ENDIF
ENDDO
ENDDO
ENDDO
@@ -134,10 +136,11 @@ SUBROUTINE compute_nadiab_moments
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qe_taue, qi_taui
+ USE model, ONLY : qe_taue, qi_taui, KIN_E
implicit none
! Add non-adiabatique term
+ IF(KIN_E) THEN
DO ip=ipsg_e,ipeg_e
IF(parray_e(ip) .EQ. 0) THEN
DO ij=ijsg_e,ijeg_e
@@ -150,6 +153,7 @@ SUBROUTINE compute_nadiab_moments
= moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
ENDIF
ENDDO
+ ENDIF
! Add non-adiabatique term
DO ip=ipsg_i,ipeg_i
IF(parray_i(ip) .EQ. 0) THEN
@@ -171,20 +175,22 @@ SUBROUTINE compute_density
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i
+ USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
IMPLICIT NONE
- IF( (ips_i .EQ. 1) .AND. (ips_e .EQ. 1) ) THEN
+ IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
DO iz = izs,ize
+ IF(KIN_E) THEN
! electron density
dens_e(ikx,iky,iz) = 0._dp
DO ij = ijs_e, ije_e
dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky,iz) * &
(moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz))
ENDDO
+ ENDIF
! ion density
dens_i(ikx,iky,iz) = 0._dp
DO ij = ijs_i, ije_i
@@ -195,6 +201,7 @@ SUBROUTINE compute_density
ENDDO
ENDDO
ENDIF
+ IF(KIN_E)&
CALL manual_3D_bcast(dens_e(ikxs:ikxe,ikys:ikye,izs:ize))
CALL manual_3D_bcast(dens_i(ikxs:ikxe,ikys:ikye,izs:ize))
END SUBROUTINE compute_density
@@ -204,7 +211,7 @@ SUBROUTINE compute_temperature
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : temp_e, temp_i, kernel_e, kernel_i
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i
+ USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
IMPLICIT NONE
REAL(dp) :: j_dp
COMPLEX(dp) :: Tperp, Tpar
@@ -215,6 +222,7 @@ SUBROUTINE compute_temperature
DO ikx = ikxs,ikxe
DO iz = izs,ize
! electron temperature
+ IF(KIN_E) THEN
temp_e(ikx,iky,iz) = 0._dp
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
@@ -223,7 +231,7 @@ SUBROUTINE compute_temperature
* (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+ SQRT2/3._dp * kernel_e(ij,ikx,iky,iz) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel)
ENDDO
-
+ ENDIF
! ion temperature
temp_i(ikx,iky,iz) = 0._dp
DO ij = ijs_i, ije_i
@@ -237,6 +245,7 @@ SUBROUTINE compute_temperature
ENDDO
ENDDO
ENDIF
+ IF(KIN_E)&
CALL manual_3D_bcast(temp_e(ikxs:ikxe,ikys:ikye,izs:ize))
CALL manual_3D_bcast(temp_i(ikxs:ikxe,ikys:ikye,izs:ize))
END SUBROUTINE compute_temperature
diff --git a/src/srcinfo.h b/src/srcinfo.h
index be0f13ec768c2ec4de3124f8bd979d5d0d90988b..85bcb0f4f0277904b7ad93de01a72c066912e435 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='59f9a2d-dirty')
+parameter (VERSION='d29af00-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Wed Oct 27 15:24:35 CEST 2021')
+parameter (EXECDATE='Fri Oct 29 18:03:02 CEST 2021')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index be0f13ec768c2ec4de3124f8bd979d5d0d90988b..85bcb0f4f0277904b7ad93de01a72c066912e435 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='59f9a2d-dirty')
+parameter (VERSION='d29af00-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Wed Oct 27 15:24:35 CEST 2021')
+parameter (EXECDATE='Fri Oct 29 18:03:02 CEST 2021')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index 26b3df805abeca601bc861086f29fd7cb3953ace..a5b65bb5efa875033fd36651faf561057e5cba21 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -9,7 +9,7 @@ SUBROUTINE stepon
USE initial_par, ONLY: WIPE_ZF, WIPE_TURB
USE ghosts
USE grid
- USE model
+ USE model, ONLY : NON_LIN, KIN_E
use prec_const
USE time_integration
USE numerics, ONLY: wipe_zonalflow, wipe_turbulence
@@ -25,7 +25,7 @@ SUBROUTINE stepon
!----- BEFORE: All fields are updated for step = n
! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
- CALL moments_eq_rhs_e
+ 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)
@@ -69,11 +69,13 @@ SUBROUTINE stepon
mlend=.FALSE.
IF(.NOT.nlend) THEN
mlend=mlend .or. checkfield(phi,' phi')
+ IF(KIN_E) THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
mlend=mlend .or. checkfield(moments_e(ip,ij,:,:,:,updatetlevel),' moments_e')
ENDDO
ENDDO
+ ENDIF
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
mlend=mlend .or. checkfield(moments_i(ip,ij,:,:,:,updatetlevel),' moments_i')
@@ -88,6 +90,7 @@ SUBROUTINE stepon
END SUBROUTINE checkfield_all
SUBROUTINE anti_aliasing
+ IF(KIN_E)THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO ikx=ikxs,ikxe
@@ -99,6 +102,7 @@ SUBROUTINE stepon
END DO
END DO
END DO
+ ENDIF
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
DO ikx=ikxs,ikxe
@@ -115,6 +119,7 @@ SUBROUTINE stepon
SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
IF ( contains_kx0 ) THEN
! Electron moments
+ IF(KIN_E) THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO iz=izs,ize
@@ -126,6 +131,7 @@ SUBROUTINE stepon
END DO
END DO
END DO
+ ENDIF
! Ion moments
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index a295575e5464fdf0efdf415e2f9f08f84fd7dc0c..556ed5c21220a822da087c70aa4bed14d051f5f9 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -4,7 +4,8 @@ MODULE utility
use prec_const
IMPLICIT NONE
- PUBLIC :: manual_2D_bcast, manual_3D_bcast
+ PUBLIC :: manual_2D_bcast, manual_3D_bcast,&
+ simpson_rule_z, o2e_z, e2o_z
CONTAINS
@@ -149,4 +150,74 @@ SUBROUTINE manual_3D_bcast(field_)
ENDIF
END SUBROUTINE manual_3D_bcast
+SUBROUTINE simpson_rule_z(f,intf)
+ ! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
+ !from molix BJ Frei
+ use prec_const
+ use grid
+ !
+ implicit none
+ !
+ complex(dp),dimension(izs:ize), intent(in) :: f
+ COMPLEX(dp), intent(out) :: intf
+ !
+ COMPLEX(dp) :: buff_
+ !
+ IF(Nz .GT. 1) THEN
+ IF(mod(Nz,2) .ne. 0 ) THEN
+ ERROR STOP 'Simpson rule: Nz must be an even number !!!!'
+ ENDIF
+ !
+ buff_ = 0._dp
+ !
+ DO iz = izs, Nz/2
+ IF(iz .eq. Nz/2) THEN ! ... iz = ize
+ buff_ = buff_ + (f(izs) + 4._dp*f(ize) + f(ize-1 ))
+ ELSE
+ buff_ = buff_ + (f(2*iz+1) + 4._dp*f(2*iz) + f(2*iz-1 ))
+ ENDIF
+ ENDDO
+ !
+ !
+ intf = buff_*deltaz/3._dp
+ !
+ ELSE
+ intf = f(izs)
+ ENDIF
+END SUBROUTINE simpson_rule_z
+
+SUBROUTINE o2e_z(fo,fe)
+ ! gives the value of a field from the odd grid to the even one
+ use prec_const
+ use grid
+ !
+ implicit none
+ !
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: fe !
+ !
+ DO iz = 2,Nz
+ fe(iz) = 0.5_dp*(fo(iz)+fo(iz-1))
+ ENDDO
+ ! Periodic boundary conditions
+ fe(1) = 0.5_dp*(fo(1) + fo(Nz))
+END SUBROUTINE o2e_z
+
+SUBROUTINE e2o_z(fe,fo)
+ ! gives the value of a field from the even grid to the odd one
+ use prec_const
+ use grid
+ !
+ implicit none
+ !
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: fo
+ !
+ DO iz = 1,Nz-1
+ fo(iz) = 0.5_dp*(fe(iz+1)+fe(iz))
+ ENDDO
+ ! Periodic boundary conditions
+ fo(Nz) = 0.5_dp*(fe(1) + fe(Nz))
+END SUBROUTINE e2o_z
+
END MODULE utility
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index cf0111b17a3296d86ad545a3480f962f04ed9258..eff2e01673390d6acd0e33b22cbfbb27d21e645c 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -4,7 +4,10 @@ outfile ='';
%% Directory of the simulation
if 1% Local results
outfile ='';
-outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
+outfile ='';
+outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_lin';
+% outfile ='fluxtube_salphaB_s0/64x64x16_5x3_L_200_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK';
+% outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
% outfile ='Linear_Device/64x64x20_5x2_Lx_20_Ly_150_q0_25_kN_0.24_kT_0.03_nu_1e-02_DGGK';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
@@ -12,7 +15,7 @@ outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
@@ -21,7 +24,7 @@ end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 03; JOBNUMMAX = 03;
+JOBNUMMIN = 00; JOBNUMMAX = 20;
% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A
compile_results %Compile the results from first output found to JOBNUMMAX if existing
@@ -33,42 +36,17 @@ default_plots_options
disp('Plots')
FMT = '.fig';
-if 0
-%% Time evolutions and growth rate
-plot_time_evolution_and_gr
-end
-
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
TAVG_0 = 1500; TAVG_1 = 4000; % Averaging times duration
plot_radial_transport_and_shear
end
-if 0
-%% Space time diagramms
-cmax = 0.00001 % max of the colorbar for transport
-tstart = 0; tend = Ts3D(end); % time window
-plot_space_time_diagrams
-end
-
-if 0
-%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 1000; tend = 4000;
-% tstart = 10000; tend = 12000;
-% Chose the field to plot
-% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
-% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0;
-plot_kperp_spectrum
-end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-t0 =000; iz = 1; ix = 1; iy = 1;
+t0 =000; iz = 2; ix = 1; iy = 1;
skip_ =1; FPS = 30; DELAY = 1/FPS;
[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
@@ -80,15 +58,15 @@ INTERP = 0; NORMALIZED = 1; CONST_CMAP = 0; BWR =1;% Gif options
% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
-FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-% FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
+FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
% FIELD = phi-Z_phi; NAME = 'NZphi'; FIELDNAME = '\phi_{NZ}';
% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
% Sliced
-% plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
+plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z';
% plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z';
-plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
+% plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y';
% K-space
% FIELD = PHI; NAME = 'PHI'; FIELDNAME = '\tilde \phi';
@@ -112,11 +90,11 @@ create_gif
end
if 0
-%% Photomaton : real space
+%% 2D plot : real space
% Chose the field to plot
% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}';
-FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
+% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i';
% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e';
% FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}';
@@ -124,16 +102,17 @@ FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}'
% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
-% FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
+FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
% FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}';
% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
% Chose when to plot it
-tf = 1500:500:3000;
+tf = [0 15 27 28 30];
-% Sliced
-ix = 1; iy = 1; iz = 1;
+% Planar plot: choose a plane to plot at x0/y0/z0 coordinates
+x0 = 0.0; y0 = 0.0; z0 = 0.0;
+[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0));
% plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')']
% plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')']
plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(z=',num2str(round(z(iz)/pi)),'\pi)']
@@ -162,35 +141,44 @@ save_figure
end
if 0
-%% Photomaton : k space
+%% 2D plot : k space
% Chose the field to plot
% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
% Chose when to plot it
-tf = 1500:500:3000;
+tf = [0 15 27 28 30];
% tf = 8000;
-% Sliced
-ix = 1; iy = 1; iz = 1;
+% Planar plot: choose a plane to plot at x0/y0/z0 coordinates
+x0 = 0.0; y0 = 0.3; z0 = 0.5*pi;
+[~,ix] = min(abs(x-x0)); [~,iy] = min(abs(y-y0)); [~,iz] = min(abs(z-z0));
% plt = @(x,it) abs(x(ix, :, :,it)); X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_x=',num2str(round(kx(ix))),')'];
% plt = @(x,it) abs(x( :,iy, :,it)); X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_y=',num2str(round(ky(iy))),')'];
-plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)'];
+% plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)'];
+% %
+% plt_x = @(x) fftshift(x,1); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,1); plt = @(x,it) max(abs(x( :, :,:,it)),[],1);
+% X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_x)'];
+%
+% plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,1); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],2);
+% X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(\max_y)'];
+
+plt_x = @(x) fftshift(x,2); plt_y = @(x) fftshift(x,2); plt_z = @(x) fftshift(x,2); plt = @(x,it) max(abs(x( :, :,:,it)),[],3);
+X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(\max_z)'];
%
TNAME = [];
fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 300*numel(tf), 500])
-plt_2 = @(x) (fftshift(x,2));
for i_ = 1:numel(tf)
[~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
subplot(1,numel(tf),i_)
DATA = plt_2(squeeze(plt(FIELD,it)));
- pclr = pcolor(fftshift(X,2),fftshift(Y,2),DATA); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
- caxis([0 1]*5e3);
+ pclr = pcolor(plt_x(X),plt_y(Y),plt_z(DATA)); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
+% caxis([0 1]*5e3);
% caxis([-1 1]*5);
xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
if(i_ > 1); set(gca,'ytick',[]); end;
@@ -227,7 +215,7 @@ xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax);
xlabel('$x/\rho_s$'); ylabel('$s_{E\times B,x}$'); grid on
end
-if 1
+if 0
%% zonal vs nonzonal energies for phi(t)
it0 = 01; itend = Ns3D;
trange = it0:itend;
@@ -310,4 +298,54 @@ plt_2 = @(x) (fftshift(x,2));
end
title(['$',FIELDLTX,'$ Zonal Spectrum']); legend('show');
save_figure
+end
+
+if 0
+%% Time evolutions and growth rate
+plot_time_evolution_and_gr
+end
+
+if 0
+%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
+% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
+tstart = 1000; tend = 4000;
+% tstart = 10000; tend = 12000;
+% Chose the field to plot
+% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
+% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
+LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0;
+plot_kperp_spectrum
+end
+
+if 0
+%% Torus plot
+aminor = EPS; % Torus minor radius
+Rmajor = 1.; % Torus major radius
+theta = linspace(-pi, pi, 30) ; % Poloidal angle
+phi = linspace(0., 2.*pi, 30) ; % Toroidal angle
+[t, p] = meshgrid(phi, theta);
+x = (Rmajor + aminor.*cos(p)) .* cos(t);
+y = (Rmajor + aminor.*cos(p)) .* sin(t);
+z = aminor.*sin(p);
+figure;
+torus=surf(x, y, z); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
+set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none')
+xlabel('X');ylabel('Y');zlabel('Z');
+% field line plot
+Nturns = 1;
+theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
+xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0);
+yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0);
+zFL = aminor.*sin(theta);
+plot3(xFL,yFL,zFL,'r')
+% Planes plot
+theta = linspace(-pi, pi, Nz) ; % Poloidal angle
+xFL = (Rmajor + aminor.*cos(theta)) .* cos(theta*Q0);
+yFL = (Rmajor + aminor.*cos(theta)) .* sin(theta*Q0);
+zFL = aminor.*sin(theta);
+plot3(xFL,yFL,zFL,'sb')
+axis equal
+
end
\ No newline at end of file
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index e1cd8d836911bdd83e512743e6e74bc04a3c0d08..95365b4d7f0b5d162b03eeef37cd6a8802cd5d5f 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -13,14 +13,14 @@ K_T = 0.0; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-Nx = 100; % real space x-gridpoints
-Ny = 1; % '' y-gridpoints
-Lx = 150; % Size of the squared frequency domain
-Ly = 1; % Size of the squared frequency domain
-Nz = 1; % number of perpendicular planes (parallel grid)
-q0 = 1.0; % safety factor
-shear = 0.0; % magnetic shear
-eps = 0.0; % inverse aspect ratio
+NX = 100; % real space x-gridpoints
+NY = 1; % '' y-gridpoints
+LX = 150; % Size of the squared frequency domain
+LY = 1; % Size of the squared frequency domain
+NZ = 1; % number of perpendicular planes (parallel grid)
+Q0 = 1.0; % safety factor
+SHEAR = 0.0; % magnetic shear
+EPS = 0.0; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 100; % Maximal time unit
DT = 1e-2; % Time step
@@ -31,8 +31,9 @@ SPS5D = 1; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-SIMID = 'test_4.1'; % Name of the simulation
+SIMID = 'Linear_entropy_mode'; % Name of the simulation
NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+KIN_E = 1;
% Collision operator
% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
CO = 2;
@@ -51,7 +52,7 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-kmax = Nx*pi/Lx;% Highest fourier mode
+kmax = NX*pi/LX;% Highest fourier mode
NU_HYP = 0.0; % Hyperdiffusivity coefficient
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
@@ -76,9 +77,9 @@ mup_ = MU_P;
muj_ = MU_J;
Nparam = numel(PA);
param_name = 'PJ';
-gamma_Ni00 = zeros(Nparam,floor(Nx/2)+1);
-gamma_Nipj = zeros(Nparam,floor(Nx/2)+1);
-gamma_phi = zeros(Nparam,floor(Nx/2)+1);
+gamma_Ni00 = zeros(Nparam,floor(NX/2)+1);
+gamma_Nipj = zeros(Nparam,floor(NX/2)+1);
+gamma_phi = zeros(Nparam,floor(NX/2)+1);
for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
@@ -96,7 +97,7 @@ for i = 1:Nparam
%%
filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
load_results
- for ikx = 1:Nx/2+1
+ for ikx = 1:NX/2+1
tend = max(Ts3D(abs(Ni00(ikx,1,1,:))~=0));
tstart = 0.6*tend;
[~,itstart] = min(abs(Ts3D-tstart));
@@ -138,10 +139,9 @@ plt = @(x) x;
end
grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(kx)]);
title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
-% title(['$\nabla N = 0$', ', $\nu=',num2str(NU),'$'])
legend('show'); %xlim([0.01,10])
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']);
+saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']);
end
end
if 0
diff --git a/wk/local_run.m b/wk/local_run.m
index 0120ca24cf0adabcc1671a10acc075916c2d65bf..0ed013d84270b3b9311354d6fc521ad08bddf5d8 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -5,30 +5,32 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
-K_N = 1/0.6; % Density gradient drive (R/Ln)
-K_T = 0.00; % Temperature gradient
-K_E = 0.00; % Electrostat gradient
+K_N = 2.22; % Density gradient drive
+K_T = 6.0; % Temperature '''
+K_E = 0.00; % Electrostat gradient
+SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NU_HYP = 0.0;
+KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
-Nx = 1024; % Spatial radial resolution ( = 2x radial modes)
-Lx = 120; % Radial window size
-Ny = 256; % Spatial azimuthal resolution (= azim modes)
-Ly = 120; % Azimuthal window size
-Nz = 1; % number of perpendicular planes (parallel grid)
-q0 = 1.0; % safety factor (Lz = 2*pi*q0)
-P = 2;
-J = 1;
+NX = 50; % Spatial radial resolution ( = 2x radial modes)
+LX = 300; % Radial window size
+NY = 100; % Spatial azimuthal resolution (= azim modes)
+LY = 300; % Azimuthal window size
+NZ = 20; % number of perpendicular planes (parallel grid)
+P = 4;
+J = 2;
%% GEOMETRY PARAMETERS
-shear = 0.0; % magnetic shear
-eps = 0.0; % inverse aspect ratio (controls parallel magnetic gradient)
-gradB = 0.0; % Magnetic gradient
-curvB = 0.0; % Magnetic curvature
+Q0 = 2.7; % safety factor
+SHEAR = 0.0; % magnetic shear
+EPS = 0.18; % inverse aspect ratio
+GRADB = 1.0; % Magnetic gradient
+CURVB = 1.0; % Magnetic curvature
%% TIME PARAMETERS
-TMAX = 90; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 10; % Maximal time unit
+DT = 5e-3; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 1/2; % Sampling per time unit for 3D arrays
+SPS3D = 5; % Sampling per time unit for 3D arrays
SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
JOB2LOAD= -1;
@@ -38,12 +40,13 @@ JOB2LOAD= -1;
CO = 1;
CLOS = 0; % Closure model (0: =0 truncation)
NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID = 'Linear_Device'; % Name of the simulation
-SIMID = 'simulation_A'; % Name of the simulation
+SIMID = 'fluxtube_salphaB_s0'; % Name of the simulation
+% SIMID = 'simulation_A'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
-NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
+NON_LIN = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
-INIT_ZF = 0; ZF_AMP = 0.0;
+INIT_PHI= 1; % Start simulation with a noisy phi (0= noisy moments 00)
+INIT_ZF = 0; ZF_AMP = 0.0;
INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
%% OUTPUTS
W_DOUBLE = 1;
@@ -62,13 +65,12 @@ KERN = 0; % Kernel model (0 : GK)
KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = Nx*pi/L;% Highest fourier mode
+kmax = NX*pi/LX;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
BCKGD0 = 0.0; % Init background
TAU = 1.0; % e/i temperature ratio
-INIT_PHI= 1; % Start simulation with a noisy phi and moments
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% Setup and file management
diff --git a/wk/plot_phi_ballooning.m b/wk/plot_phi_ballooning.m
index 47bd05e8dd54710f288ef043796e77075d89e0d7..678cd3b9f891c3f53a5d9978a160ecc2ac3dc3db 100644
--- a/wk/plot_phi_ballooning.m
+++ b/wk/plot_phi_ballooning.m
@@ -4,14 +4,14 @@ phi_imag=(imag(PHI(:,:,:,it)));
% Apply baollooning tranform
for iky=2
dims = size(phi_real);
- phib_real = zeros( dims(1)*Nz ,1);
+ phib_real = zeros( dims(1)*dims(3) ,1);
phib_imag= phib_real;
b_angle = phib_real;
- midpoint = floor((dims(1)*Nz )/2)+1;
+ midpoint = floor((dims(1)*dims(3) )/2)+1;
for ip =1: dims(1)
- start_ = (ip -1)*Nz +1;
- end_ = ip*Nz;
+ start_ = (ip -1)*dims(3) +1;
+ end_ = ip*dims(3);
phib_real(start_:end_) = phi_real(ip,iky,:);
phib_imag(start_:end_) = phi_imag(ip,iky, :);
end
@@ -20,8 +20,8 @@ for iky=2
Nkx = numel(kx)-1; coordz = z;
idx = -Nkx:1:Nkx;
for ip =1: dims(1)
- for iz=1:Nz
- ii = Nz*(ip -1) + iz;
+ for iz=1:dims(3)
+ ii = dims(3)*(ip -1) + iz;
b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
end
end
@@ -30,9 +30,9 @@ for iky=2
[~,idxLFS] = min(abs(b_angle -0));
phib = phib_real(:) + 1i * phib_imag(:);
% Normalize to the outermid plane
- phib_norm = phib(:);% / phib( idxLFS) ;
- phib_real_norm(:) = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
- phib_imag_norm(:) = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
+ phib_norm = phib;% / phib( idxLFS) ;
+ phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
+ phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
figure; hold all;
@@ -43,26 +43,4 @@ for iky=2
xlabel('$\chi / \pi$')
ylabel('$\phi_B (\chi)$');
title(['HeLaZ(-) molix(o) benchmark, t=',num2str(Ts3D(it))]);
-% title(['HeLaZ,$(P,J) =(',num2str(PMAXI),', ', num2str(JMAXI),'$), $\nu =',num2str(NU),...
-% '$, $\epsilon = ',num2str(eps),'$, $k_y = ', num2str(ky( iky)),'$, $q =',num2str(q0),'$, $s = ', num2str(shear),'$, $R_N = ', ...
-% num2str(K_N),'$, $R_{T_i} = ', num2str(K_T),'$, $N_z =',num2str(Nz),'$']);
- %set(gca,'Yscale','log')
- %
-
-% %Check symmetry of the mode at the outter mid plane
-% figure; hold all;
-% right = phib_real(midpoint+1:end);
-% left = fliplr(phib_real(1:midpoint-1)');
-% up_down_symm = right(1:end) - left(1:end-1)';
-% %plot(b_angle(midpoint+1:end)/pi,phib_real(midpoint+1:end),'-xb');
-% plot(b_angle(midpoint+1:end)/pi,up_down_symm ,'-xb');
- %plot(abs(b_angle(1:midpoint-1))/pi,phib_real(1:midpoint-1),'-xb');
- %
- %
- % figure; hold all
- % plot(b_angle/pi, phib_imag.^2 + phib_real.^2 ,'xk');
- % %set(gca,'Yscale','log')
- % xlabel('$\chi / \pi$')
- % ylabel('$\phi_B (\chi)$');
- % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),'$, $\epsilon = ',num2str(epsilon),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $k_y =',num2str(ky),'$']);
end
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
index 4fa96d088e66daa73d613648ac11afefab9f8931..8700441432f73a1f6e4f8678eaa887ffc7a7184f 100644
--- a/wk/shearless_linear_fluxtube.m
+++ b/wk/shearless_linear_fluxtube.m
@@ -11,15 +11,16 @@ TAU = 1.0; % e/i temperature ratio
K_N = 2.22; % Density gradient drive
K_T = 6.0; % Temperature '''
SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 1; % Kinetic (1) or adiabatic (0) electron model
%% GRID PARAMETERS
-Nx = 1; % real space x-gridpoints
-Ny = 2; % '' y-gridpoints
-Lx = 0; % Size of the squared frequency domain
-Ly = 2*pi/0.25; % Size of the squared frequency domain
-Nz = 24; % number of perpendicular planes (parallel grid)
-q0 = 2.7; % safety factor
-shear = 0.0; % magnetic shear
-eps = 0.18; % inverse aspect ratio
+NX = 1; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
+LX = 0; % Size of the squared frequency domain
+LY = 2*pi/0.25; % Size of the squared frequency domain
+NZ = 24; % number of perpendicular planes (parallel grid)
+Q0 = 2.7; % safety factor
+SHEAR = 0.0; % magnetic shear
+EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 10; % Maximal time unit
DT = 1e-3; % Time step
@@ -90,7 +91,7 @@ for i = 1:Nparam
% system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0 > out.txt; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
% Load and process results
%%