diff --git a/src/cosolver_interface_mod.F90 b/src/cosolver_interface_mod.F90
index a96bceaf804cf47d47b88fceb92a29f952107251..888c722ef140c6ccd6f2cfd7eedb097ccf9a7ad3 100644
--- a/src/cosolver_interface_mod.F90
+++ b/src/cosolver_interface_mod.F90
@@ -13,9 +13,8 @@ CONTAINS
USE parallel, ONLY: num_procs_p, comm_p,dsp_p,rcv_p
USE grid, ONLY: &
local_na, &
- local_np, total_np, parray, ngp,&
- total_nj,jarray, ngj,&
- local_nkx, local_nky, local_nz, ngz, bar
+ local_np, total_np, total_nj,&
+ local_nkx, local_nky, local_nz, bar
USE array, ONLY: Capj
USE MPI
IMPLICIT NONE
@@ -24,19 +23,13 @@ CONTAINS
COMPLEX(dp), DIMENSION(local_np) :: TColl_distr
COMPLEX(dp) :: Tmp_
INTEGER :: iz,ikx,iky,ij,ip,ia,ikx_C,iky_C,iz_C
- INTEGER :: p_int,j_int, ierr
- INTEGER :: ipi, iji, izi
+ INTEGER :: ierr
z:DO iz = 1,local_nz
- izi = iz+ngz/2
x:DO ikx = 1,local_nkx
y:DO iky = 1,local_nky
a:DO ia = 1,local_na
j:DO ij = 1,total_nj
- iji = ij+ngj/2
- j_int = jarray(iji)
p:DO ip = 1,total_np
- ipi = ip + ngp/2
- p_int = parray(ipi)
!! Take GK or DK limit
IF (GK_CO) THEN ! GK operator (k-dependant)
ikx_C = ikx; iky_C = iky; iz_C = iz;
@@ -75,8 +68,8 @@ CONTAINS
SUBROUTINE apply_cosolver_mat(ia,ip,ij,iky,ikx,iz,ikx_C,iky_C,iz_C,local_sum)
USE grid, ONLY: &
local_na, &
- local_np, parray, ngp,&
- total_nj, jarray, dmax, ngj, bar, ngz
+ local_np, parray,parray_full, ngp,&
+ total_nj, jarray,jarray_full, dmax, ngj, bar, ngz
USE array, ONLY: Caa, Cab_T, Cab_F, nadiab_moments
USE model, ONLY: CLOS
USE species, ONLY: nu_ab
@@ -85,12 +78,10 @@ CONTAINS
COMPLEX(dp), INTENT(OUT) :: local_sum
INTEGER :: ib,iq,il
INTEGER :: p_int,q_int,j_int,l_int
- INTEGER :: ipi, iji, izi, iqi, ili
+ INTEGER :: izi, iqi, ili
izi = iz+ngz/2
- iji = ij+ngj/2
- j_int = jarray(iji)
- ipi = ip + ngp/2
- p_int = parray(ipi)
+ p_int = parray_full(ip)
+ j_int = jarray_full(ij)
!! Apply the cosolver collision matrix
local_sum = 0._dp ! Initialization
q:DO iq = 1,local_np
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index ab4c2515bbecad9048a352939fdebecd507ef7a3..c9861d1c4788cb6eb756a3a8359282b501676a32 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -63,16 +63,16 @@ SUBROUTINE update_ghosts_p_mom
USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
ngp,ngj,ngz
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(local_np+ngp) :: slice_z
+ COMPLEX(dp), DIMENSION(local_np+ngp) :: slice_p
INTEGER :: ia,iz,ij,iy,ix
- DO iz = 1,local_nz
+ DO iz = 1+ngz/2,local_nz+ngz/2
DO ix = 1,local_nkx
DO iy = 1,local_nky
- DO ij = 1,local_nj
+ DO ij = 1+ngj/2,local_nj+ngj/2
DO ia = 1,local_na
- slice_z = moments(ia,:,ij,iy,ix,iz,updatetlevel)
- CALL exchange_ghosts_1D(slice_z,nbr_L,nbr_R,local_np,ngp)
- moments(ia,:,ij,iy,ix,iz,updatetlevel) = slice_z
+ slice_p = moments(ia,:,ij,iy,ix,iz,updatetlevel)
+ CALL exchange_ghosts_1D(slice_p,nbr_L,nbr_R,local_np,ngp)
+ moments(ia,:,ij,iy,ix,iz,updatetlevel) = slice_p
ENDDO
ENDDO
ENDDO
@@ -210,7 +210,7 @@ SUBROUTINE update_ghosts_z_3D(field)
DO iky = 1,local_nky
DO ikx = 1,local_nkx
ikxBC_L = ikx_zBC_L(iky,ikx);
- IF (ikxBC_L .GT. 0) THEN ! Exchanging the modes that have a periodic pair (a)
+ IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
DO ig = 1,ngz/2
field(iky,ikx,first-ig) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-ig)
ENDDO
@@ -220,7 +220,7 @@ SUBROUTINE update_ghosts_z_3D(field)
ENDDO
ENDIF
ikxBC_R = ikx_zBC_R(iky,ikx);
- IF (ikxBC_R .GT. 0) THEN ! Exchanging the modes that have a periodic pair (a)
+ IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
DO ig = 1,ngz/2
field(iky,ikx,last+ig) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,ig)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index f24ed835477543ec4b48467f19d787a3bf4a46e9..74d1615b26e922471e225ce385481f095445d5fc 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -123,15 +123,15 @@ SUBROUTINE compute_moments_eq_rhs
RHS = RHS - mu_p*diff_p_coeff*p_int**6*Napj
IF (j_int .GT. 1) &
RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
- ! CASE('dvpar4')
- ! ! fourth order numerical diffusion in vpar
- ! IF(p_int .GE. 4) &
- ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! ! (not used often so not parallelized)
- ! RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
- ! ! + dummy Laguerre diff
- ! IF (j_int .GE. 2) &
- ! RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
+ CASE('dvpar4')
+ ! fourth order numerical diffusion in vpar
+ IF(p_int .GE. 4) &
+ ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! (not used often so not parallelized)
+ RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
+ ! + dummy Laguerre diff
+ IF (j_int .GE. 2) &
+ RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
CASE DEFAULT
END SELECT
ELSE
@@ -163,8 +163,6 @@ SUBROUTINE compute_moments_eq_rhs
print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
print*,'---'
IF(updatetlevel .EQ. 4) stop
- ! print*, xnapjp1(2,:)
- ! print*, pp2
stop
END SUBROUTINE compute_moments_eq_rhs
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index e8f3aac55d2e94731f4e5e5d10e83fdba5de7af1..e37be6cad785028aee2aebddf34240d4f1785b40 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -137,9 +137,9 @@ SUBROUTINE evaluate_poisson_op
USE model, ONLY : ADIAB_E, tau_e
USE prec_const, ONLY: dp
IMPLICIT NONE
- REAL(dp) :: pol_ion, pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
+ REAL(dp) :: pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
INTEGER :: in,ikx,iky,iz,ia
- REAL(dp) :: pol_i, pol_e, sumker, spol ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
+ REAL(dp) :: sumker ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
! This term has no staggered grid dependence. It is evalued for the
! even z grid since poisson uses p=0 moments and phi only.
@@ -180,37 +180,37 @@ END SUBROUTINE evaluate_poisson_op
SUBROUTINE evaluate_ampere_op
USE prec_const, ONLY : dp
USE array, ONLY : kernel, inv_ampere_op
- USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, ngz, &
- jmax, kparray, kxarray, kyarray, SOLVE_AMPERE, iodd
+ USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, ngz, total_nj, ngj,&
+ kparray, kxarray, kyarray, SOLVE_AMPERE, iodd
USE model, ONLY : beta
USE species, ONLY : q, sigma
USE geometry, ONLY : hatB
USE prec_const, ONLY: dp
IMPLICIT NONE
- REAL(dp) :: pol_tot, kperp2 ! (Z^2/tau (1-sum_n kernel_na^2))
+ REAL(dp) :: sum_jpol, kperp2, operator ! (Z^2/tau (1-sum_n kernel_na^2))
INTEGER :: in,ikx,iky,iz,ia
! We do not solve Ampere if beta = 0 to spare waste of ressources
IF(SOLVE_AMPERE) THEN
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO iz = 1,local_nz
- kperp2 = kparray(iky,ikx,iz,1)**2
+ x:DO ikx = 1,local_nkx
+ y:DO iky = 1,local_nky
+ z:DO iz = 1,local_nz
+ kperp2 = kparray(iky,ikx,iz+ngz/2,iodd)**2
IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
inv_ampere_op(iky, ikx, iz) = 0._dp
ELSE
- !!!!!!!!!!!!!!!!! Ion contribution
- pol_tot = 0._dp
- DO ia = 1,local_na
+ sum_jpol = 0._dp
+ a:DO ia = 1,local_na
! loop over n only up to the max polynomial degree
- DO in=1,jmax+1
- pol_tot = pol_tot + q(ia)**2/(sigma(ia)**2)*kernel(ia,in,iky,ikx,iz,1)**2 ! ... sum recursively ...
- END DO
- END DO
- inv_ampere_op(iky, ikx, iz) = 1._dp/(2._dp*kperp2*hatB(iz+ngz/2,iodd)**2 + beta*pol_tot)
+ j:DO in=1,total_nj
+ sum_jpol = sum_jpol + q(ia)**2/(sigma(ia)**2)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,iodd)**2 ! ... sum recursively ...
+ END DO j
+ END DO a
+ operator = 2._dp*kperp2*hatB(iz+ngz/2,iodd)**2 + beta*sum_jpol
+ inv_ampere_op(iky, ikx, iz) = 1._dp/operator
ENDIF
- END DO
- END DO
- END DO
+ END DO z
+ END DO y
+ END DO x
ENDIF
END SUBROUTINE evaluate_ampere_op
!******************************************************************************!
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 5cc934350aed6dbc2bee9adfb18b5c5e48e731e9..b38d9cd13ef5d4ed18e9e73e64c09c95c917c488 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -84,10 +84,9 @@ CONTAINS
call simpson_rule_z(local_nz,deltaz,integrant,integral)
! Get process local gyrocenter flux with a factor two to account for the negative ky modes
gflux_local = 2._dp*integral*iInt_Jacobian
-
!
- integrant = 0._dp ! reset auxiliary variable
!!---------- Particle flux (gyrokinetic) ------------
+ integrant = 0._dp ! reset auxiliary variable
! Electrostatic part
IF(CONTAINSp0) THEN
DO iz = 1,local_nz ! we take interior points only
@@ -121,14 +120,9 @@ CONTAINS
ENDDO
ENDIF
! Integrate over z
- ! print*, integrant
call simpson_rule_z(local_nz,deltaz,integrant,integral)
- ! print*, integral
- ! stop
! Get process local particle flux with a factor two to account for the negative ky modes
pflux_local = 2._dp*integral*iInt_Jacobian
- print*,REAL(pflux_local,dp)
-
!!!!---------- Sum over all processes ----------
buffer(1) = REAL(gflux_local,dp)
buffer(2) = REAL(pflux_local,dp)
@@ -196,7 +190,7 @@ CONTAINS
DO iky = 1,local_nky
DO in = 1, local_nj
ini = in + ngj/2 !interior index for ghosted arrays
- n_dp = jarray(in)
+ n_dp = jarray(ini)
integrant(iz) = integrant(iz) &
+Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
*kernel(ia,ini,iky,ikx,izi,iodd)*(&
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 8f2120a3b0f2518c3c489230e6142e74c20e22d1..db3b3dd453c5bd3fe6479d2a90a3c44833df16eb 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -107,25 +107,27 @@ CONTAINS
use species, ONLY: sqrt_tau_o_sigma, q
use model, ONLY: beta
IMPLICIT NONE
- COMPLEX(dp) :: iota ! current density
- INTEGER :: in, ia, ikx, iky, iz
+ COMPLEX(dp) :: j_a ! current density
+ INTEGER :: in, ia, ikx, iky, iz, ini, izi
! Execution time start
CALL cpu_time(t0_poisson)
!! Ampere can be solved only with beta > 0 and for process containng p=1 moments
IF ( SOLVE_AMPERE ) THEN
z:DO iz = 1,local_nz
- x:DO ikx = 1,local_nky
- y:DO iky = 1,local_nkx
- !!!!!!!!!!!!!!! compute current density contribution "iota = q_a u_a" for each species
- iota = 0._dp
+ izi = iz+ngz/2
+ x:DO ikx = 1,local_nkx
+ y:DO iky = 1,local_nky
+ !!!!!!!!!!!!!!! compute current density contribution "j_a = q_a u_a" for each species
+ j_a = 0._dp
n:DO in=1,total_nj
+ ini = in+ngj/2
a:DO ia = 1,local_na
- iota = iota &
- +q(ia)*sqrt_tau_o_sigma(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,in,iky,ikx,iz+ngz/2,updatetlevel)
+ j_a = j_a &
+ +q(ia)*sqrt_tau_o_sigma(ia)*kernel(ia,ini,iky,ikx,izi,iodd)*moments(ia,ip1,ini,iky,ikx,izi,updatetlevel)
ENDDO a
ENDDO n
!!!!!!!!!!!!!!! Inverting the Ampere equation
- psi(iky,ikx,iz+ngz/2) = beta*inv_ampere_op(iky,ikx,iz)*iota
+ psi(iky,ikx,iz+ngz/2) = beta*inv_ampere_op(iky,ikx,iz)*j_a
ENDDO y
ENDDO x
ENDDO z
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index 9437edb8162bb303f50c00acf4eb2ef4f8459119..f6a8807c823adec78718240e0ca647988b504acd 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -1,5 +1,5 @@
&BASIC
- nrun = 1
+ nrun = 99999999
dt = 0.01
tmax = 5
maxruntime = 356400
@@ -8,9 +8,9 @@
&GRID
pmax = 4
jmax = 1
- Nx = 4
+ Nx = 2
Lx = 200
- Ny = 2
+ Ny = 4
Ly = 60
Nz = 6
SG = .f.
@@ -31,11 +31,11 @@
shift_y= 0.0
/
&OUTPUT_PAR
- dtsave_0d = 0.01
+ dtsave_0d = 0.5
dtsave_1d = -1
dtsave_2d = -1
- dtsave_3d = 0.01
- dtsave_5d = 0.01
+ dtsave_3d = 1
+ dtsave_5d = 5
write_doubleprecision = .f.
write_gamma = .t.
write_hf = .t.
@@ -58,9 +58,9 @@
mu_z = 0.6
HYP_V = 'hypcoll'
mu_p = 0.1
- mu_j = 0.
+ mu_j = 0.5
nu = 1.0
- beta = 0
+ beta = 0.1
ADIAB_E = .f.
tau_e = 1.0
/
@@ -77,14 +77,14 @@
! electrons
name_ = 'electrons'
tau_ = 1.0
- sigma_= 1!0.023338
- q_ = 1!-1.0
+ sigma_= 0.023338
+ q_ = -1.0
k_N_ = 1.0!2.22
k_T_ = 2.0!6.96
/
&COLLISION_PAR
- collision_model = 'SG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index 1b515d515953d085d3ba50570066b8d215dccd06..d4f127f005d537334dcb0124ec23614f0fc78dba 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -1,5 +1,5 @@
&BASIC
- nrun = 1
+ nrun = 99999999
dt = 0.01
tmax = 5
maxruntime = 356400
@@ -9,9 +9,9 @@
jmaxe = 1
pmaxi = 4
jmaxi = 1
- Nx = 4
+ Nx = 2
Lx = 200
- Ny = 2
+ Ny = 4
Ly = 60
Nz = 6
SG = .f.
@@ -24,11 +24,11 @@
parallel_bc = 'dirichlet'
/
&OUTPUT_PAR
- nsave_0d = 1
+ nsave_0d = 50
nsave_1d = -1
nsave_2d = -1
- nsave_3d = 1
- nsave_5d = 1
+ nsave_3d = 100
+ nsave_5d = 500
write_doubleprecision = .f.
write_gamma = .t.
write_hf = .t.
@@ -51,22 +51,22 @@
mu_z = 0.6
HYP_V = 'hypcoll'
mu_p = 0.1
- mu_j = 0.
+ mu_j = 0.5
nu = 1.0
tau_e = 1
tau_i = 1
- sigma_e = 1!0.023338
+ sigma_e = 0.023338
sigma_i = 1
- q_e = 1!-1
+ q_e = -1
q_i = 1
K_Ne = 1!6.96
K_Te = 2!2.22
K_Ni = 3!6.96
K_Ti = 4!2.22
- beta = 0
+ beta = 0.1
/
&COLLISION_PAR
- collision_model = 'SG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
gyrokin_CO = .t.
interspecies = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'