diff --git a/Makefile b/Makefile
index 80ae86717d3d9a6419f5f2b8f0560935a45693f9..9112278d95986de1ce17654219df6c1fc245ed03 100644
--- a/Makefile
+++ b/Makefile
@@ -61,16 +61,18 @@ mvmod:
$(OBJDIR)/diagnose.o : src/srcinfo.h
-FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
-$(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/collision_mod.o \
-$(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
-$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
-$(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/geometry_mod.o \
-$(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs_mod.o \
-$(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o \
-$(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o \
-$(OBJDIR)/restarts_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o \
-$(OBJDIR)/utility_mod.o
+FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o \
+$(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o \
+$(OBJDIR)/collision_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o \
+$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \
+$(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/geometry_mod.o \
+$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/inital.o \
+$(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
+$(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o \
+$(OBJDIR)/parallel_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o \
+$(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o \
+$(OBJDIR)/restarts_mod.o $(OBJDIR)/solve_EM_fields.o $(OBJDIR)/stepon.o \
+$(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(EXEC): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
@@ -82,7 +84,7 @@ $(OBJDIR)/utility_mod.o
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
$(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o \
- $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o \
$(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
@@ -114,11 +116,6 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
- $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o\
- $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
-
$(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/geometry_mod.o \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o \
$(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
@@ -149,7 +146,7 @@ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
$(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
$(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o\
@@ -157,7 +154,9 @@ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
$(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o \
+ $(OBJDIR)/solve_EM_fields.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
$(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o \
@@ -176,18 +175,19 @@ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
$(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 $(OBJDIR)/array_mod.o \
- $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@
+ $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
+ $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o\
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
+
$(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
$(OBJDIR)/coeff_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics_mod.F90 -o $@
- $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o \
- $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/parallel_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
-
$(OBJDIR)/parallel_mod.o : src/parallel_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
$(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/parallel_mod.F90 -o $@
@@ -213,11 +213,19 @@ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
$(OBJDIR)/restarts_mod.o : src/restarts_mod.F90 $(OBJDIR)/diagnostics_par_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/restarts_mod.F90 -o $@
+ $(OBJDIR)/solve_EM_fields.o : src/solve_EM_fields.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o \
+ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/parallel_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/solve_EM_fields.F90 -o $@
+
$(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+ $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o \
+ $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
$(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
diff --git a/matlab/setup.m b/matlab/setup.m
index 6288c1de688fdc68d47fe40347b7544c096bf008..6228cac6a33725440c569c406db5704100c8259f 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -23,6 +23,7 @@ MODEL.NL_CLOS = NL_CLOS;
MODEL.LINEARITY = ['''',LINEARITY,''''];
MODEL.KIN_E = KIN_E;
if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
+MODEL.beta = BETA;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
MODEL.N_HD = N_HD;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 17838bc320df017124f422c0407357e41b41e1ee..90dee735fcacc54f79262fe2fa2ba42e544ac7e3 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -71,9 +71,10 @@ fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
fprintf(fid,[' K_n = ', num2str(MODEL.K_n),'\n']);
fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']);
fprintf(fid,[' K_E = ', num2str(MODEL.K_E),'\n']);
-fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
-fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
+fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
+fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
+fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&COLLISION_PAR\n');
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 2db05d5edc043bc8991184ea3efe888441e2368a..8cb3a2ef147732d93cfd495002d289b6a952464d 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -56,13 +56,15 @@ MODULE array
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1 ! mirror lin coeff for adiab mom
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_e, xphijp1_e, xphijm1_e
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_i, xphijp1_i, xphijm1_i
-
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xpsij_e, xpsijp1_e, xpsijm1_e
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xpsij_i, xpsijp1_i, xpsijm1_i
! Kernel function evaluation (ij,ikx,iky,iz,odd/even p)
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
! Poisson operator (ikx,iky,iz)
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_ampere_op
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 50a0042b12b43576c955909f53380b71696f4e7e..5df1f8c087a13cc3c65500220c13282eb0ec7790 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -209,6 +209,7 @@ SUBROUTINE diagnose_full(kstep)
CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
+ IF (write_phi) CALL creatg(fidres, "/data/var3d/psi", "psi")
IF (write_Na00) THEN
IF(KIN_E)&
@@ -405,6 +406,7 @@ SUBROUTINE diagnose_3d
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
+ IF (write_phi) CALL write_field3d_kykxz(psi (ikys:ikye,ikxs:ikxe,izs:ize), 'psi')
IF (write_Na00) THEN
IF(KIN_E)THEN
diff --git a/src/fields_mod.F90 b/src/fields_mod.F90
index eac8c04f227a1e86df6cb05b938339c0a7dafc17..728b9fdbf47df65964cb6b509c63300689f125c6 100644
--- a/src/fields_mod.F90
+++ b/src/fields_mod.F90
@@ -12,4 +12,7 @@ MODULE fields
! Normalized electric potential: \hat{\phi} ! (kx,ky,z)
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: phi
+!------------------Vector field part
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: psi
+
END MODULE fields
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 2fa4bc65dba3ab3b3d88fa159f6f80f34222d42c..1121680a848598cc73c941edb4ec290fcd931bf7 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -1,15 +1,14 @@
module ghosts
USE basic
-USE fields, ONLY : moments_e, moments_i, phi
USE grid
USE time_integration
-USE model, ONLY : KIN_E
+USE model, ONLY : KIN_E, beta
USE geometry, ONLY : SHEARED, ikx_zBC_L, ikx_zBC_R
IMPLICIT NONE
INTEGER :: status(MPI_STATUS_SIZE), source, dest, count, ipg
-PUBLIC :: update_ghosts_moments, update_ghosts_phi
+PUBLIC :: update_ghosts_moments, update_ghosts_EM
CONTAINS
@@ -31,15 +30,19 @@ SUBROUTINE update_ghosts_moments
tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
END SUBROUTINE update_ghosts_moments
-SUBROUTINE update_ghosts_phi
+SUBROUTINE update_ghosts_EM
CALL cpu_time(t0_ghost)
IF(Nz .GT. 1) THEN
CALL update_ghosts_z_phi
+
+ IF(beta .GT. 0._dp) &
+ CALL update_ghosts_z_psi
ENDIF
tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
-END SUBROUTINE update_ghosts_phi
+END SUBROUTINE update_ghosts_EM
+
!Communicate p+1, p+2 moments to left neighboor and p-1, p-2 moments to right one
! [a b|C D|e f] : proc n has moments a to f where a,b,e,f are ghosts
@@ -54,7 +57,7 @@ END SUBROUTINE update_ghosts_phi
! ^ ^
!Closure by zero truncation : 0 0
SUBROUTINE update_ghosts_p_e
-
+ USE fields, ONLY : moments_e
IMPLICIT NONE
count = (ijge_e-ijgs_e+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1)
@@ -82,7 +85,7 @@ END SUBROUTINE update_ghosts_p_e
!Communicate p+1, p+2 moments to left neighboor and p-1, p-2 moments to right one
SUBROUTINE update_ghosts_p_i
-
+ USE fields, ONLY : moments_i
IMPLICIT NONE
count = (ijge_i-ijgs_i+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1) ! Number of elements sent
@@ -121,6 +124,7 @@ END SUBROUTINE update_ghosts_p_i
!Periodic: 0 1
SUBROUTINE update_ghosts_z_e
USE parallel, ONLY : buff_pjxy_zBC_e
+ USE fields, ONLY : moments_e
IMPLICIT NONE
INTEGER :: ikxBC_L, ikxBC_R
IF(Nz .GT. 1) THEN
@@ -178,6 +182,7 @@ END SUBROUTINE update_ghosts_z_e
SUBROUTINE update_ghosts_z_i
USE parallel, ONLY : buff_pjxy_zBC_i
+ USE fields, ONLY : moments_i
IMPLICIT NONE
INTEGER :: ikxBC_L, ikxBC_R
IF(Nz .GT. 1) THEN
@@ -235,6 +240,7 @@ END SUBROUTINE update_ghosts_z_i
SUBROUTINE update_ghosts_z_phi
USE parallel, ONLY : buff_xy_zBC
+ USE fields, ONLY : phi
IMPLICIT NONE
INTEGER :: ikxBC_L, ikxBC_R
CALL cpu_time(t1_ghost)
@@ -294,5 +300,67 @@ SUBROUTINE update_ghosts_z_phi
tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
END SUBROUTINE update_ghosts_z_phi
+SUBROUTINE update_ghosts_z_psi
+ USE parallel, ONLY : buff_xy_zBC
+ USE fields, ONLY : psi
+ IMPLICIT NONE
+ INTEGER :: ikxBC_L, ikxBC_R
+ CALL cpu_time(t1_ghost)
+ IF(Nz .GT. 1) THEN
+ IF (num_procs_z .GT. 1) THEN
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ count = (ikye-ikys+1) * (ikxe-ikxs+1)
+ !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv( psi(:,:,ize ), count, MPI_DOUBLE_COMPLEX, nbr_U, 30, & ! Send to Up the last
+ buff_xy_zBC(:,:,-1), count, MPI_DOUBLE_COMPLEX, nbr_D, 30, & ! Receive from Down the first-1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv( psi(:,:,ize-1), count, MPI_DOUBLE_COMPLEX, nbr_U, 31, & ! Send to Up the last
+ buff_xy_zBC(:,:,-2), count, MPI_DOUBLE_COMPLEX, nbr_D, 31, & ! Receive from Down the first-2
+ comm0, status, ierr)
+
+ !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv( psi(:,:,izs ), count, MPI_DOUBLE_COMPLEX, nbr_D, 32, & ! Send to Down the first
+ buff_xy_zBC(:,:,+1), count, MPI_DOUBLE_COMPLEX, nbr_U, 32, & ! Recieve from Up the last+1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv( psi(:,:,izs+1), count, MPI_DOUBLE_COMPLEX, nbr_D, 33, & ! Send to Down the first
+ buff_xy_zBC(:,:,+2), count, MPI_DOUBLE_COMPLEX, nbr_U, 33, & ! Recieve from Up the last+2
+ comm0, status, ierr)
+ ELSE
+ buff_xy_zBC(:,:,-1) = psi(:,:,ize )
+ buff_xy_zBC(:,:,-2) = psi(:,:,ize-1)
+ buff_xy_zBC(:,:,+1) = psi(:,:,izs )
+ buff_xy_zBC(:,:,+2) = psi(:,:,izs+1)
+ ENDIF
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ ikxBC_L = ikx_zBC_L(iky,ikx);
+ IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
+ ! first-1 gets last
+ psi(iky,ikx,izs-1) = buff_xy_zBC(iky,ikxBC_L,-1)
+ ! first-2 gets last-1
+ psi(iky,ikx,izs-2) = buff_xy_zBC(iky,ikxBC_L,-2)
+ ELSE
+ psi(iky,ikx,izs-1) = 0._dp
+ psi(iky,ikx,izs-2) = 0._dp
+ ENDIF
+ ikxBC_R = ikx_zBC_R(iky,ikx);
+ IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
+ ! last+1 gets first
+ psi(iky,ikx,ize+1) = buff_xy_zBC(iky,ikxBC_R,+1)
+ ! last+2 gets first+1
+ psi(iky,ikx,ize+2) = buff_xy_zBC(iky,ikxBC_R,+2)
+ ELSE
+ psi(iky,ikx,ize+1) = 0._dp
+ psi(iky,ikx,ize+2) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+ CALL cpu_time(t1_ghost)
+ tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
+END SUBROUTINE update_ghosts_z_psi
+
END MODULE ghosts
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 2120042ff98ee876c620c840692ee206c1b7445e..a31128d017ac0141b858a8cfd4d47fe05657c722 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -99,6 +99,7 @@ MODULE grid
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip0_e, CONTAINS_ip0_i
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip1_e, CONTAINS_ip1_i
LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip2_e, CONTAINS_ip2_i
+ LOGICAL, PUBLIC, PROTECTED :: SOLVE_POISSON, SOLVE_AMPERE
! Usefull inverse numbers
REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz
@@ -125,25 +126,14 @@ CONTAINS
Nx, Lx, Ny, Ly, Nz, Npol, SG
READ(lu_in,grid)
+ IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use
+ SG = .FALSE.
+
!! Compute the maximal degree of full GF moments set
! i.e. : all moments N_a^pj s.t. p+2j<=d are simulated (see GF closure)
dmaxe = min(pmaxe,2*jmaxe+1)
dmaxi = min(pmaxi,2*jmaxi+1)
- ! If no parallel dim (Nz=1), the moment hierarchy is separable between odds and even P
- !! and since the energy is injected in P=0 and P=2 for density/temperature gradients
- !! there is no need of simulating the odd p which will only be damped.
- !! We define in this case a grid Parray = 0,2,4,...,Pmax i.e. deltap = 2 instead of 1
- !! to spare computation
- IF(Nz .EQ. 1) THEN
- deltape = 2; deltapi = 2;
- pp2 = 1; ! index p+2 is ip+1
- SG = .FALSE.
- ELSE
- deltape = 1; deltapi = 1;
- pp2 = 2; ! index p+2 is ip+1
- ENDIF
-
! Usefull precomputations
inv_Nx = 1._dp/REAL(Nx,dp)
inv_Ny = 1._dp/REAL(Ny,dp)
@@ -160,9 +150,24 @@ CONTAINS
SUBROUTINE set_pgrid
USE prec_const
+ USE model, ONLY: beta ! To know if we solve ampere or not and put odd p moments
IMPLICIT NONE
INTEGER :: ip, istart, iend, in
+ ! If no parallel dim (Nz=1) and no EM effects (beta=0), the moment hierarchy
+ !! is separable between odds and even P and since the energy is injected in
+ !! P=0 and P=2 for density/temperature gradients there is no need of
+ !! simulating the odd p which will only be damped.
+ !! We define in this case a grid Parray = 0,2,4,...,Pmax i.e. deltap = 2
+ !! instead of 1 to spare computation
+ IF((Nz .EQ. 1) .AND. (beta .EQ. 0._dp)) THEN
+ deltape = 2; deltapi = 2;
+ pp2 = 1; ! index p+2 is ip+1
+ ELSE
+ deltape = 1; deltapi = 1;
+ pp2 = 2; ! index p+2 is ip+2
+ ENDIF
+
! Total number of Hermite polynomials we will evolve
total_np_e = (Pmaxe/deltape) + 1
total_np_i = (Pmaxi/deltapi) + 1
@@ -204,6 +209,8 @@ CONTAINS
CONTAINS_ip0_i = .FALSE.
CONTAINS_ip1_i = .FALSE.
CONTAINS_ip2_i = .FALSE.
+ SOLVE_POISSON = .FALSE.
+ SOLVE_AMPERE = .FALSE.
ALLOCATE(parray_e(ipgs_e:ipge_e))
ALLOCATE(parray_i(ipgs_i:ipge_i))
DO ip = ipgs_e,ipge_e
@@ -212,10 +219,12 @@ CONTAINS
IF(parray_e(ip) .EQ. 0) THEN
ip0_e = ip
CONTAINS_ip0_e = .TRUE.
+ SOLVE_POISSON = .TRUE.
ENDIF
IF(parray_e(ip) .EQ. 1) THEN
ip1_e = ip
CONTAINS_ip1_e = .TRUE.
+ SOLVE_AMPERE = .TRUE.
ENDIF
IF(parray_e(ip) .EQ. 2) THEN
ip2_e = ip
@@ -228,10 +237,12 @@ CONTAINS
IF(parray_i(ip) .EQ. 0) THEN
ip0_i = ip
CONTAINS_ip0_i = .TRUE.
- ENDIF
+ SOLVE_POISSON = .TRUE.
+ ENDIF
IF(parray_i(ip) .EQ. 1) THEN
ip1_i = ip
CONTAINS_ip1_i = .TRUE.
+ SOLVE_AMPERE = .TRUE.
ENDIF
IF(parray_i(ip) .EQ. 2) THEN
ip2_i = ip
@@ -245,6 +256,11 @@ CONTAINS
! Precomputations
pmaxe_dp = real(pmaxe,dp)
pmaxi_dp = real(pmaxi,dp)
+
+ ! Overwrite SOLVE_AMPERE flag if beta is zero
+ IF(beta .EQ. 0._dp) THEN
+ SOLVE_AMPERE = .FALSE.
+ ENDIF
END SUBROUTINE set_pgrid
SUBROUTINE set_jgrid
diff --git a/src/inital.F90 b/src/inital.F90
index c692fb10565f77aef3b3c5f589b1b91985208179..d1958e72f83b388214b488aa40e17c9cd9bdac39 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -8,7 +8,7 @@ SUBROUTINE inital
USE time_integration, ONLY: set_updatetlevel
USE collision, ONLY: load_COSOlver_mat, cosolver_coll
USE closure, ONLY: apply_closure_model
- USE ghosts, ONLY: update_ghosts_moments, update_ghosts_phi
+ USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
USE restarts, ONLY: load_moments, job2load
USE numerics, ONLY: play_with_modes, save_EM_ZF_modes
USE processing, ONLY: compute_fluid_moments
@@ -24,8 +24,8 @@ SUBROUTINE inital
IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
CALL load_moments ! get N_0
CALL update_ghosts_moments
- CALL poisson ! compute phi_0=phi(N_0)
- CALL update_ghosts_phi
+ CALL solve_EM_fields ! compute phi_0=phi(N_0)
+ CALL update_ghosts_EM
! through initialization
ELSE
SELECT CASE (INIT_OPT)
@@ -33,35 +33,35 @@ SUBROUTINE inital
CASE ('phi')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
- CALL update_ghosts_phi
+ CALL update_ghosts_EM
! set moments_00 (GC density) with noise and compute phi afterwards
CASE('mom00')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
CALL init_gyrodens ! init only gyrocenter density
CALL update_ghosts_moments
- CALL poisson
- CALL update_ghosts_phi
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
! init all moments randomly (unadvised)
CASE('allmom')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments'
CALL init_moments ! init all moments
CALL update_ghosts_moments
- CALL poisson
- CALL update_ghosts_phi
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
! init a gaussian blob in gyrodens
CASE('blob')
IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
CALL initialize_blob
CALL update_ghosts_moments
- CALL poisson
- CALL update_ghosts_phi
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
! init moments 00 with a power law similarly to GENE
CASE('ppj')
IF (my_id .EQ. 0) WRITE(*,*) 'ppj init ~ GENE'
call init_ppj
CALL update_ghosts_moments
- CALL poisson
- CALL update_ghosts_phi
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
END SELECT
ENDIF
! closure of j>J, p>P and j<0, p<0 moments
@@ -70,7 +70,7 @@ SUBROUTINE inital
! ghosts for p parallelization
IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
CALL update_ghosts_moments
- CALL update_ghosts_phi
+ CALL update_ghosts_EM
!! End of phi and moments initialization
! Save (kx,0) and (0,ky) modes for num exp
@@ -399,8 +399,8 @@ SUBROUTINE init_ppj
REAL(dp) :: kx, ky, sigma_z, amp, ky_shift, z
INTEGER, DIMENSION(12) :: iseedarr
- sigma_z = pi/4.0
- amp = 1e4
+ sigma_z = pi/4._dp
+ amp = 1.0_dp
!**** Broad noise initialization *******************************************
! Electrons
@@ -418,11 +418,11 @@ SUBROUTINE init_ppj
IF(ky .EQ. 0) THEN
moments_e(ip,ij,iky,ikx,iz,:) = 0._dp
ELSE
- moments_e(ip,ij,iky,ikx,iz,:) = 0._dp!0.5_dp * ky_min/(ABS(ky)+ky_min)
+ moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
IF(ky .GT. 0) THEN
- moments_e(ip,ij,iky,ikx,iz,:) = 0._dp!(kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ moments_e(ip,ij,iky,ikx,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(kx_min/(ABS(kx)+kx_min))
ENDIF
@@ -461,11 +461,11 @@ SUBROUTINE init_ppj
IF(ky .EQ. 0) THEN
moments_i(ip,ij,iky,ikx,iz,:) = 0._dp
ELSE
- moments_i(ip,ij,iky,ikx,iz,:) = 0._dp!0.5_dp * ky_min/(ABS(ky)+ky_min)
+ moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
IF(ky .GT. 0) THEN
- moments_i(ip,ij,iky,ikx,iz,:) = 0._dp!(kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ moments_i(ip,ij,iky,ikx,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(kx_min/(ABS(kx)+kx_min))
ENDIF
@@ -509,9 +509,5 @@ SUBROUTINE init_ppj
ENDDO
ENDIF
- ! Adjust the scaling to trigger faster NL saturation
- IF(KIN_E) &
- moments_e = 1e3*moments_e
- moments_i = 1e3*moments_i
END SUBROUTINE init_ppj
!******************************************************************************!
diff --git a/src/memory.F90 b/src/memory.F90
index f39ca92dc316eda624558556b5e62782981f5b44..ec0297e87d7ae04fced4ed2d0725093a9288fa73 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -14,9 +14,11 @@ SUBROUTINE memory
! Electrostatic potential
CALL allocate_array( phi, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( psi, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
CALL allocate_array( phi_EM, ikys,ikye, izs,ize)
CALL allocate_array(inv_poisson_op, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( inv_ampere_op, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array( inv_pol_ion, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array(HF_phi_correction_operator, ikys,ikye, ikxs,ikxe, izs,ize)
@@ -100,10 +102,16 @@ SUBROUTINE memory
CALL allocate_array( xphij_e, ips_e,ipe_e, ijs_e,ije_e)
CALL allocate_array( xphijp1_e, ips_e,ipe_e, ijs_e,ije_e)
CALL allocate_array( xphijm1_e, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( xpsij_e, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( xpsijp1_e, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( xpsijm1_e, ips_e,ipe_e, ijs_e,ije_e)
ENDIF
CALL allocate_array( xphij_i, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijp1_i, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijm1_i, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( xpsij_i, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( xpsijp1_i, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( xpsijm1_i, ips_i,ipe_i, ijs_i,ije_i)
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 043c5b02da8ce832b9e45381678865884a8d453f..12d16a584c117061a45101ae009e2246dca5e291 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -30,6 +30,7 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient
REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature
REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length
+ REAL(dp), PUBLIC, PROTECTED :: beta = 0._dp ! electron plasma Beta (8piNT_e/B0^2)
REAL(dp), PUBLIC, PROTECTED :: taue_qe ! factor of the magnetic moment coupling
REAL(dp), PUBLIC, PROTECTED :: taui_qi !
@@ -46,24 +47,30 @@ MODULE model
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
-
+ REAL(dp), PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i
+ REAL(dp), PUBLIC, PROTECTED :: sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i
PUBLIC :: model_readinputs, model_outputinputs
CONTAINS
SUBROUTINE model_readinputs
! Read the input parameters
- USE basic, ONLY : lu_in
+ USE basic, ONLY : lu_in, my_id
USE prec_const
IMPLICIT NONE
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
mu_x, mu_y, N_HD, mu_z, 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
+ K_n, K_T, K_E, GradB, CurvB, lambdaD, beta
READ(lu_in,model_par)
+ IF(.NOT. KIN_E) THEN
+ IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0'
+ beta = 0._dp
+ ENDIF
+
taue_qe = tau_e/q_e ! factor of the magnetic moment coupling
taui_qi = tau_i/q_i ! factor of the magnetic moment coupling
qe_taue = q_e/tau_e
@@ -78,6 +85,10 @@ CONTAINS
sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
sqrt_sigmae2_taue_o2 = SQRT(sigma_e**2 * tau_e/2._dp) ! to avoid multiple SQRT eval
sqrt_sigmai2_taui_o2 = SQRT(sigma_i**2 * tau_i/2._dp)
+ q_o_sqrt_tau_sigma_e = q_e/SQRT(tau_e)/sigma_e ! For psi field terms
+ q_o_sqrt_tau_sigma_i = q_i/SQRT(tau_i)/sigma_i ! For psi field terms
+ sqrt_tau_o_sigma_e = SQRT(tau_e)/sigma_e ! For Ampere eq
+ sqrt_tau_o_sigma_i = SQRT(tau_i)/sigma_i
!! We use the ion-ion collision as normalization with definition
! nu_ii = 4 sqrt(pi)/3 T_i^(-3/2) m_i^(-1/2) q^4 n_i0 ln(Lambda)
!
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index d6ce46df38dfe838f26b04eb363efa3e0cb5c184..9d1626b254f8e40532cefce100422eb8f540b32c 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -35,9 +35,9 @@ SUBROUTINE moments_eq_rhs_e
REAL(dp) :: kx, ky, kperp2, dzlnB_o_J
COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives
COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
- COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
+ COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi
COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: i_kx,i_ky,phikykxz
+ COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -52,6 +52,7 @@ SUBROUTINE moments_eq_rhs_e
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
phikykxz = phi(iky,ikx,iz)! tmp phi value
+ psikykxz = psi(iky,ikx,iz)! tmp psi value
! Kinetic loops
jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1
@@ -93,13 +94,22 @@ SUBROUTINE moments_eq_rhs_e
Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = (xphij_i (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
- + xphijp1_i(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
- + xphijm1_i(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz
+ Tphi = (xphij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
+ + xphijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
+ + xphijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz
ELSE
Tphi = 0._dp
ENDIF
+ !! Vector potential term
+ IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3
+ Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
+ + xpsijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
+ + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz
+ ELSE
+ Tpsi = 0._dp
+ ENDIF
+
!! Sum of all RHS terms
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
@@ -109,7 +119,7 @@ SUBROUTINE moments_eq_rhs_e
! Mirror term (parallel magnetic gradient)
- gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) &
! Drives (density + temperature gradients)
- - i_ky * Tphi &
+ - i_ky * (Tphi - Tpsi) &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
- mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
@@ -165,8 +175,8 @@ SUBROUTINE moments_eq_rhs_i
COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1
COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
- COMPLEX(dp) :: i_kx, i_ky, phikykxz
+ COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi
+ COMPLEX(dp) :: i_kx, i_ky, phikykxz, psikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -181,6 +191,7 @@ SUBROUTINE moments_eq_rhs_i
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
phikykxz = phi(iky,ikx,iz)! tmp phi value
+ psikykxz = psi(iky,ikx,iz)! tmp phi value
! Kinetic loops
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
@@ -231,6 +242,16 @@ SUBROUTINE moments_eq_rhs_i
Tphi = 0._dp
ENDIF
+ !! Vector potential term
+ IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3
+ Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) &
+ + xpsijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) &
+ + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz
+ ELSE
+ Tpsi = 0._dp
+ ENDIF
+
+
!! Sum of all RHS terms
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
@@ -240,7 +261,7 @@ SUBROUTINE moments_eq_rhs_i
! Mirror term (parallel magnetic gradient)
- gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir &
! Drives (density + temperature gradients)
- - i_ky * Tphi &
+ - i_ky * (Tphi - Tpsi) &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
- mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 833c4f680f6fb1a147c49c108aff04638b82c197..1ede5ac6f9b392e781d09ea75deb627fb7891e06 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -6,8 +6,8 @@ MODULE numerics
USE coeff
implicit none
- PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, compute_lin_coeff
- PUBLIC :: play_with_modes, save_EM_ZF_modes
+ PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, evaluate_ampere_op
+ PUBLIC :: compute_lin_coeff, play_with_modes, save_EM_ZF_modes
CONTAINS
@@ -151,10 +151,57 @@ SUBROUTINE evaluate_poisson_op
END SUBROUTINE evaluate_poisson_op
!******************************************************************************!
+!******************************************************************************!
+!!!!!!! Evaluate inverse polarisation operator for Poisson equation
+!******************************************************************************!
+SUBROUTINE evaluate_ampere_op
+ USE basic
+ USE array, Only : kernel_e, kernel_i, inv_ampere_op
+ USE grid
+ USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E, beta
+ IMPLICIT NONE
+ REAL(dp) :: pol_i, pol_e, kperp2 ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
+ INTEGER :: ini,ine
+
+ ! We do not solve Ampere if beta = 0 to spare waste of ressources
+ IF(SOLVE_AMPERE) THEN
+ ! This term has no staggered grid dependence. It is evalued for the
+ ! even z grid since poisson uses p=0 moments and phi only.
+ kxloop: DO ikx = ikxs,ikxe
+ kyloop: DO iky = ikys,ikye
+ zloop: DO iz = izs,ize
+ kperp2 = kparray(iky,ikx,iz,0)**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
+ ! loop over n only if the max polynomial degree
+ pol_i = 0._dp
+ DO ini=1,jmaxi+1
+ pol_i = pol_i + kernel_i(ini,iky,ikx,iz,0)**2 ! ... sum recursively ...
+ END DO
+ pol_i = q_i**2/(sigma_i**2) * pol_i
+ !!!!!!!!!!!!! Electron contribution
+ pol_e = 0._dp
+ ! loop over n only if the max polynomial degree
+ DO ine=1,jmaxe+1 ! ine = n+1
+ pol_e = pol_e + kernel_e(ine,iky,ikx,iz,0)**2 ! ... sum recursively ...
+ END DO
+ pol_e = q_e**2/(sigma_e**2) * pol_e
+ inv_ampere_op(iky, ikx, iz) = 1._dp/(2._dp*kperp2 + beta*(pol_i + pol_e))
+ ENDIF
+ END DO zloop
+ END DO kyloop
+ END DO kxloop
+ ENDIF
+
+END SUBROUTINE evaluate_ampere_op
+!******************************************************************************!
+
SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
- K_T, K_n, CurvB, GradB, KIN_E
+ K_T, K_n, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i
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
@@ -270,7 +317,7 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_i(ip,ij) =+K_n + 2.*j_dp*K_T
+ xphij_i(ip,ij) =+K_n + 2._dp*j_dp*K_T
xphijp1_i(ip,ij) =-K_T*(j_dp+1._dp)
xphijm1_i(ip,ij) =-K_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
@@ -282,6 +329,48 @@ SUBROUTINE compute_lin_coeff
ENDIF
ENDDO
ENDDO
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !! EM linear coefficients for moment RHS !!!!!!!!!!
+ IF (KIN_E) THEN
+ DO ip = ips_e, ipe_e
+ p_int= parray_e(ip) ! Hermite degree
+ DO ij = ijs_e, ije_e
+ j_int= jarray_e(ij) ! REALof Laguerre degree
+ j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+ !! Electrostatic potential pj terms
+ IF (p_int .EQ. 1) THEN ! kronecker p1
+ xpsij_e(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_e)/sigma_e
+ xpsijp1_e(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
+ xpsijm1_e(ip,ij) =- K_T* j_dp * SQRT(tau_e)/sigma_e
+ ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
+ xpsij_e(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
+ xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp;
+ ELSE
+ xpsij_e(ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp
+ xpsijm1_e(ip,ij) = 0._dp;
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+ DO ip = ips_i, ipe_i
+ p_int= parray_i(ip) ! Hermite degree
+ DO ij = ijs_i, ije_i
+ j_int= jarray_i(ij) ! REALof Laguerre degree
+ j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+ !! Electrostatic potential pj terms
+ IF (p_int .EQ. 1) THEN ! kronecker p1
+ xpsij_i(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_i)/sigma_i
+ xpsijp1_i(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
+ xpsijm1_i(ip,ij) =- K_T* j_dp * SQRT(tau_i)/sigma_i
+ ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
+ xpsij_i(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
+ xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp;
+ ELSE
+ xpsij_i(ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp
+ xpsijm1_i(ip,ij) = 0._dp;
+ ENDIF
+ ENDDO
+ ENDDO
END SUBROUTINE compute_lin_coeff
!******************************************************************************!
diff --git a/src/poisson.F90 b/src/poisson.F90
deleted file mode 100644
index 8b401b0c2c053df45827e6ab847a912e96cb4f49..0000000000000000000000000000000000000000
--- a/src/poisson.F90
+++ /dev/null
@@ -1,78 +0,0 @@
-SUBROUTINE poisson
- ! Solve poisson equation to get phi
-
- USE basic
- USE time_integration, ONLY: updatetlevel
- USE array
- USE fields
- USE grid
- USE calculus, ONLY : simpson_rule_z
- USE parallel, ONLY : manual_3D_bcast
- 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
- COMPLEX(dp) :: fsa_phi, intf_ ! current flux averaged phi
- INTEGER :: count !! mpi integer to broadcast the electric potential at the end
- COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e, integrant ! charge density q_a n_a and aux var
-
- ! Execution time start
- CALL cpu_time(t0_poisson)
- !! Poisson can be solved only for process containing p=0
- IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
-
-
- kxloop: DO ikx = ikxs,ikxe
- kyloop: DO iky = ikys,ikye
- phi(iky,ikx,izs:ize) = 0._dp
- !!!! Compute ion particle charge density q_i n_i
- rho_i = 0._dp
- DO ini=1,jmaxi+1
- rho_i(izs:ize) = rho_i(izs:ize) &
- +q_i*kernel_i(ini,iky,ikx,izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)
- END DO
-
- !!!! Compute electron particle charge density q_e n_e
- rho_e = 0._dp
- IF (KIN_E) THEN ! Kinetic electrons
- DO ine=1,jmaxe+1
- rho_e(izs:ize) = rho_e(izs:ize) &
- +q_e*kernel_e(ine,iky,ikx,izs:ize,0)*moments_e(ip0_e,ine,iky,ikx,izs:ize,updatetlevel)
- END DO
- ELSE ! Adiabatic electrons
- ! Adiabatic charge density (linked to flux surface averaged phi)
- ! We compute the flux surface average solving a flux surface averaged
- ! Poisson equation, i.e.
- ! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
- ! inv_pol_ion^-1 fsa_phi = simpson(Jacobian rho_i ) * iInt_Jacobian
- fsa_phi = 0._dp
- IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
- ! Prepare integrant for z-average
- integrant(izs:ize) = Jacobian(izs:ize,0)*rho_i(izs:ize)*inv_pol_ion(iky,ikx,izs:ize)
- call simpson_rule_z(integrant(izs:ize),intf_) ! get the flux averaged phi
- fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
- ENDIF
- rho_e(izs:ize) = fsa_phi
- ENDIF
- !!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
- phi(iky,ikx,izs:ize) = (rho_e(izs:ize) + rho_i(izs:ize))*inv_poisson_op(iky,ikx,izs:ize)
-
- END DO kyloop
- END DO kxloop
-
- ! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) phi(iky_0,ikx_0,:) = 0._dp
-
- ENDIF
-
- ! Transfer phi to all the others process along p
- CALL manual_3D_bcast(phi(ikys:ikye,ikxs:ikxe,izs:ize))
-
- ! Execution time end
- CALL cpu_time(t1_poisson)
- tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-END SUBROUTINE poisson
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 649d51ef5873857790ce7e40091a37c5ca0f7425..79d80f27deeecebc1693f8147cbd3e0b839cd9fa 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -295,12 +295,13 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!
! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
!
- USE fields, ONLY : moments_i, moments_e, phi
+ USE fields, ONLY : moments_i, moments_e, phi, psi
USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, &
ddz_nepj, ddz4_Nepj, interp_nepj,&
ddz_nipj, ddz4_Nipj, interp_nipj
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qe_taue, qi_taui, KIN_E, CLOS
+ USE model, ONLY : qe_taue, qi_taui,q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i, &
+ KIN_E, CLOS
USE calculus, ONLY : grad_z, grad_z4, interp_z
IMPLICIT NONE
INTEGER :: eo, p_int, j_int
@@ -315,6 +316,9 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) &
+ qe_taue*kernel_e(ij,:,:,:,0)*phi(:,:,:)
ENDDO
+ ELSEIF(parray_e(ip) .EQ. 1) THEN
+ nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) &
+ - q_o_sqrt_tau_sigma_e*kernel_e(ij,:,:,:,0)*psi(:,:,:)
ELSE
DO ij=ijgs_e,ijge_e
nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel)
@@ -329,6 +333,9 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) &
+ qi_taui*kernel_i(ij,:,:,:,0)*phi(:,:,:)
ENDDO
+ ELSEIF(parray_i(ip) .EQ. 1) THEN
+ nadiab_moments_e(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) &
+ - q_o_sqrt_tau_sigma_i*kernel_i(ij,:,:,:,0)*psi(:,:,:)
ELSE
DO ij=ijgs_i,ijge_i
nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel)
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
new file mode 100644
index 0000000000000000000000000000000000000000..b102ac1a6fd157ccacf6918dd34d2ff6bd4f749d
--- /dev/null
+++ b/src/solve_EM_fields.F90
@@ -0,0 +1,147 @@
+SUBROUTINE solve_EM_fields
+ ! Solve Poisson and Ampere equations
+ USE model, ONLY : beta
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+
+ CALL poisson
+
+ IF(beta .GT. 0._dp) &
+ CALL ampere
+
+CONTAINS
+
+ SUBROUTINE poisson
+ ! Solve poisson equation to get phi
+ USE time_integration, ONLY: updatetlevel
+ USE array, ONLY: kernel_e, kernel_i, inv_poisson_op, inv_pol_ion
+ USE fields, ONLY: phi, moments_e, moments_i
+ USE grid
+ USE calculus, ONLY : simpson_rule_z
+ USE parallel, ONLY : manual_3D_bcast
+ use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD, KIN_E
+ USE processing, ONLY : compute_density
+ USE geometry, ONLY : iInt_Jacobian, Jacobian
+ IMPLICIT NONE
+
+ INTEGER :: ini,ine, i_, root_bcast
+ COMPLEX(dp) :: fsa_phi, intf_ ! current flux averaged phi
+ INTEGER :: count !! mpi integer to broadcast the electric potential at the end
+ COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e, integrant ! charge density q_a n_a and aux var
+
+ ! Execution time start
+ CALL cpu_time(t0_poisson)
+ !! Poisson can be solved only for process containing p=0
+ ! IF ( SOLVE_POISSON ) THEN
+ IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
+
+ kxloop: DO ikx = ikxs,ikxe
+ kyloop: DO iky = ikys,ikye
+ phi(iky,ikx,izs:ize) = 0._dp
+ !!!! Compute ion particle charge density q_i n_i
+ rho_i = 0._dp
+ DO ini=1,jmaxi+1
+ rho_i(izs:ize) = rho_i(izs:ize) &
+ +q_i*kernel_i(ini,iky,ikx,izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)
+ END DO
+
+ !!!! Compute electron particle charge density q_e n_e
+ rho_e = 0._dp
+ IF (KIN_E) THEN ! Kinetic electrons
+ DO ine=1,jmaxe+1
+ rho_e(izs:ize) = rho_e(izs:ize) &
+ +q_e*kernel_e(ine,iky,ikx,izs:ize,0)*moments_e(ip0_e,ine,iky,ikx,izs:ize,updatetlevel)
+ END DO
+ ELSE ! Adiabatic electrons
+ ! Adiabatic charge density (linked to flux surface averaged phi)
+ ! We compute the flux surface average solving a flux surface averaged
+ ! Poisson equation, i.e.
+ ! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
+ ! inv_pol_ion^-1 fsa_phi = simpson(Jacobian rho_i ) * iInt_Jacobian
+ fsa_phi = 0._dp
+ IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
+ ! Prepare integrant for z-average
+ integrant(izs:ize) = Jacobian(izs:ize,0)*rho_i(izs:ize)*inv_pol_ion(iky,ikx,izs:ize)
+ call simpson_rule_z(integrant(izs:ize),intf_) ! get the flux averaged phi
+ fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
+ ENDIF
+ rho_e(izs:ize) = fsa_phi
+ ENDIF
+ !!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
+ phi(iky,ikx,izs:ize) = (rho_e(izs:ize) + rho_i(izs:ize))*inv_poisson_op(iky,ikx,izs:ize)
+
+ END DO kyloop
+ END DO kxloop
+
+ ! Cancel origin singularity
+ IF (contains_kx0 .AND. contains_ky0) phi(iky_0,ikx_0,:) = 0._dp
+
+ ENDIF
+
+ ! Transfer phi to all the others process along p
+ CALL manual_3D_bcast(phi(ikys:ikye,ikxs:ikxe,izs:ize))
+
+ ! Execution time end
+ CALL cpu_time(t1_poisson)
+ tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ END SUBROUTINE poisson
+
+ SUBROUTINE ampere
+ ! Solve ampere equation to get psi
+ USE time_integration, ONLY: updatetlevel
+ USE array, ONLY: kernel_e, kernel_i, inv_ampere_op
+ USE fields, ONLY: moments_i, moments_e, psi
+ USE grid
+ USE parallel, ONLY : manual_3D_bcast
+ use model, ONLY : sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i, q_e, q_i, beta, KIN_E
+ IMPLICIT NONE
+
+ INTEGER :: ini,ine, i_, root_bcast
+ INTEGER :: count !! mpi integer to broadcast the electric potential at the end
+ COMPLEX(dp), DIMENSION(izs:ize) :: I_i, I_e ! current density q_a u_a and aux var
+
+ ! Execution time start
+ CALL cpu_time(t0_poisson)
+ !! Ampere can be solved only with beta > 0 and for process containing p=1
+ IF ( SOLVE_AMPERE ) THEN
+
+ kxloop: DO ikx = ikxs,ikxe
+ kyloop: DO iky = ikys,ikye
+ psi(iky,ikx,izs:ize) = 0._dp
+ !!!! Compute ion particle current density q_i u_i
+ I_i = 0._dp
+ DO ini=1,jmaxi+1
+ I_i(izs:ize) = I_i(izs:ize) &
+ +kernel_i(ini,iky,ikx,izs:ize,0)*moments_i(ip1_i,ini,iky,ikx,izs:ize,updatetlevel)
+ END DO
+ I_i = q_i * sqrt_tau_o_sigma_i * I_i
+ !!!! Compute electron particle charge density q_e n_e
+ I_e = 0._dp
+ DO ine=1,jmaxe+1
+ I_e(izs:ize) = I_e(izs:ize) &
+ +q_e*kernel_e(ine,iky,ikx,izs:ize,0)*moments_e(ip1_e,ine,iky,ikx,izs:ize,updatetlevel)
+ END DO
+ I_e = q_e * sqrt_tau_o_sigma_e * I_e
+ !!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
+ psi(iky,ikx,izs:ize) = beta*(I_e(izs:ize) + I_i(izs:ize))*inv_ampere_op(iky,ikx,izs:ize)
+
+ END DO kyloop
+ END DO kxloop
+
+ ! Cancel origin singularity
+ IF (contains_kx0 .AND. contains_ky0) psi(iky_0,ikx_0,:) = 0._dp
+
+ ENDIF
+
+ ! Transfer phi to all the others process along p
+ CALL manual_3D_bcast(psi(ikys:ikye,ikxs:ikxe,izs:ize))
+
+ ! Execution time end
+ CALL cpu_time(t1_poisson)
+ tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ END SUBROUTINE ampere
+
+END SUBROUTINE solve_EM_fields
diff --git a/src/stepon.F90 b/src/stepon.F90
index 09ad45332450af7c0b4d6d9e44c839d7e66f96c8..1771985ee4be7ff11d15e3c0a518aff934f41912 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -3,7 +3,7 @@ SUBROUTINE stepon
USE advance_field_routine, ONLY: advance_time_level, advance_field, advance_moments
USE basic
USE closure
- USE ghosts, ONLY: update_ghosts_moments, update_ghosts_phi
+ USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
USE grid
USE model, ONLY : LINEARITY, KIN_E
use prec_const
@@ -34,9 +34,9 @@ SUBROUTINE stepon
! Exchanges the ghosts values of N_n+1
CALL update_ghosts_moments
- ! Update electrostatic potential phi_n = phi(N_n+1)
- CALL poisson
- CALL update_ghosts_phi
+ ! Update electrostatic potential phi_n = phi(N_n+1) and potential vect psi
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
! Numerical experiments
! Store or cancel/maintain zonal modes artificially
@@ -130,7 +130,7 @@ SUBROUTINE stepon
END SUBROUTINE anti_aliasing
SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
- USE fields, ONLY: phi, moments_e, moments_i
+ USE fields, ONLY: phi, psi, moments_e, moments_i
IMPLICIT NONE
IF ( contains_ky0 ) THEN
! Electron moments
@@ -165,6 +165,12 @@ SUBROUTINE stepon
END DO
! must be real at origin
phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize))
+ ! Psi
+ DO ikx=2,Nkx/2 !symmetry at ky = 0
+ psi(iky_0,ikx,izs:ize) = psi(iky_0,Nkx+2-ikx,izs:ize)
+ END DO
+ ! must be real at origin
+ psi(iky_0,ikx_0,izs:ize) = REAL(psi(iky_0,ikx_0,izs:ize))
ENDIF
END SUBROUTINE enforce_symmetry
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 99ceca4dccdd71e6e1303e938173ab83738250a9..7721a364a50f299fed8bed4af9a168d08747bdac 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -44,18 +44,18 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
-options.NAME = 'v_y';
+% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'xz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
-options.COMP = 1;
+options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D(1:20:end);
+options.TIME = data.Ts3D(1:2:end);
% options.TIME = [750:1:1000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
@@ -65,7 +65,7 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = '\phi';
@@ -74,11 +74,11 @@ options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxz';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [500 800 1000];
+options.TIME = [120 150 155];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -87,9 +87,9 @@ end
if 0
%% 3D plot on the geometry
options.INTERP = 1;
-options.NAME = 'n_i';
-options.PLANES = [1:10:80];
-options.TIME = [200];
+options.NAME = '\phi';
+options.PLANES = [1:1:16];
+options.TIME = [15];
options.PLT_MTOPO = 0;
data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
@@ -121,7 +121,7 @@ options.ST = 1;
options.PLOT_TYPE = 'space-time';
options.NORMALIZED = 1;
options.JOBNUM = 0;
-options.TIME = [100:1200];
+options.TIME = [1000];
options.specie = 'i';
options.compz = 'avg';
fig = show_moments_spectrum(data,options);
@@ -167,7 +167,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
-options.TIME = data.Ts3D;
+options.TIME = [0:160];
options.NMA = 1;
options.NMODES = 15;
options.iz = 'avg';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 822daf314d50993c2084e1830c90c0b50852e4c2..4b20e0b514f73158bf4046250ca94836175bf6de 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -1,6 +1,6 @@
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
-% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.2/';
+% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/rm_corrections_HF/';
% folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
@@ -47,8 +47,8 @@ options.NAME = '\phi';
options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
-options.COMP = 'avg';
-options.TIME = [1:10];
+options.COMP = 1;
+options.TIME = [15 50 100 200];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -64,7 +64,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 7088014cefef8885ea45a1473b26152e88d1585a..9b6139f7c5a8e86ee113fe1187204fb31c77ee5d 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -16,7 +16,7 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
% outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
-%% CBC
+%% shearless CBC
% outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
% outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
% outfile ='shearless_cyclone/64x32x16x5x3_CBC_120';
@@ -31,9 +31,15 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3';
% outfile ='Zpinch_rerun/Kn_2.5_312x196x5x3_Lx_400_Ly_200';
% outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
-% outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
+outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
% outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
% outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
% outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
-outfile ='Zpinch_rerun/Kn_1.6_256x128x21x11';
+% outfile ='Zpinch_rerun/Kn_1.6_256x128x21x11';
+
+%% CBC
+% outfile = 'CBC/64x32x16x5x3';
+% outfile = 'CBC/64x128x16x5x3';
+% outfile = 'CBC/128x64x16x5x3';
+outfile = 'CBC/64x64x16x3x2';
run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 737428b07666606848abaee502fb2ca97f9184c9..33faeed8839a66ad2a6331674501d4f89b85bfc0 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -16,21 +16,22 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_N = 2.22;%2.0; % Density gradient drive
-K_T = 6.92;%0.25*K_N; % Temperature '''
+K_T = 6.96;%0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0e-1; % electron plasma beta
%% GRID PARAMETERS
-PMAXE = 4; % Hermite basis size of electrons
-JMAXE = 2; % Laguerre "
+PMAXE = 12; % Hermite basis size of electrons
+JMAXE = 6; % Laguerre "
PMAXI = 4; % " ions
JMAXI = 2; % "
-NX = 12; % real space x-gridpoints
-NY = 8; % '' y-gridpoints
+NX = 16; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.1; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
+LY = 2*pi/0.3; % Size of the squared frequency domain
+NZ = 32; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
@@ -41,7 +42,7 @@ SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 20; % Maximal time unit
-DT = 5e-3; % Time step
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
@@ -92,9 +93,9 @@ setup
% Run linear simulation
if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
@@ -107,7 +108,7 @@ JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
%% Short analysis
-if 0
+if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
nplots = 1;
@@ -118,9 +119,9 @@ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
end
-if 1
+if 0
%% Ballooning plot
-options.time_2_plot = [1]*data.Ts3D(end);
+options.time_2_plot = [120];
options.kymodes = [0.1 0.2 0.3];
options.normalized = 1;
options.field = 'phi';