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

Advances in implementation of ExB shear 

-1D routines development
-NL factor to correct moving grid method
-still an issue with the 1D r2c FFT
-ExB not working yet.
parent 9a0dd00d
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Makefile
+ 4
− 5
View file @ 1dccb6f8
...@@ -217,7 +217,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o ...@@ -217,7 +217,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
$(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 \ $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 \
$(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
$(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 \ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 \
...@@ -273,8 +273,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o ...@@ -273,8 +273,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@
$(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 \ $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 \
$(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o \ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ExB_shear_flow_mod.o \
$(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_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 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
...@@ -337,8 +337,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o ...@@ -337,8 +337,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
$(OBJDIR)/utility_mod.o : src/utility_mod.F90 \ $(OBJDIR)/utility_mod.o : src/utility_mod.F90 \
$(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
$(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
$(OBJDIR)/CLA_mod.o : src/CLA_mod.F90 \ $(OBJDIR)/CLA_mod.o : src/CLA_mod.F90 \
......
src/ExB_shear_flow_mod.F90
+ 35
− 10
View file @ 1dccb6f8
...@@ -2,15 +2,17 @@ MODULE ExB_shear_flow ...@@ -2,15 +2,17 @@ MODULE ExB_shear_flow
! This module contains the necessary tools to implement ExB shearing flow effects. ! This module contains the necessary tools to implement ExB shearing flow effects.
! The algorithm is taken from the presentation of Hammett et al. 2006 (APS) and ! The algorithm is taken from the presentation of Hammett et al. 2006 (APS) and
! it the one used in GS2. ! it the one used in GS2.
USE prec_const, ONLY: xp USE prec_const, ONLY: xp, imagu
IMPLICIT NONE IMPLICIT NONE
! Variables ! Variables
REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky) REAL(xp), PUBLIC, PROTECTED :: t0, inv_t0 = 0._xp ! charact. shear time
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jump_ExB ! jump to do to shift the kx grids REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky)
LOGICAL, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jump_ExB ! jump to do to shift the kx grids
LOGICAL, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: ExB_NL_factor! factor for nonlinear term
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: inv_ExB_NL_factor
! Routines ! Routines
PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
...@@ -18,7 +20,8 @@ CONTAINS ...@@ -18,7 +20,8 @@ CONTAINS
! Setup the variables for the ExB shear ! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow SUBROUTINE Setup_ExB_shear_flow
USE grid, ONLY : local_nky USE grid, ONLY : total_nkx, local_nky, deltakx, deltaky
USE model, ONLY : ExBrate
IMPLICIT NONE IMPLICIT NONE
! Setup the ExB shift ! Setup the ExB shift
...@@ -30,6 +33,20 @@ CONTAINS ...@@ -30,6 +33,20 @@ CONTAINS
jump_ExB = 0 jump_ExB = 0
ALLOCATE(shiftnow_ExB(local_nky)) ALLOCATE(shiftnow_ExB(local_nky))
shiftnow_ExB = .FALSE. shiftnow_ExB = .FALSE.
! Setup nonlinear factor
ALLOCATE( ExB_NL_factor(total_nkx,local_nky))
ALLOCATE(inv_ExB_NL_factor(total_nkx,local_nky))
ExB_NL_factor = 1._xp
inv_ExB_NL_factor = 1._xp
IF(ExBrate .NE. 0) THEN
t0 = deltakx/deltaky/ExBrate
inv_t0 = 1._xp/t0
ELSE ! avoid 1/0 division (t0 is killed anyway in this case)
t0 = 0._xp
inv_t0 = 0._xp
ENDIF
END SUBROUTINE Setup_ExB_shear_flow END SUBROUTINE Setup_ExB_shear_flow
! Update according to the current ExB shear value ! Update according to the current ExB shear value
...@@ -109,12 +126,14 @@ CONTAINS ...@@ -109,12 +126,14 @@ CONTAINS
! update the ExB shear value for the next time step ! update the ExB shear value for the next time step
SUBROUTINE Update_ExB_shear_flow SUBROUTINE Update_ExB_shear_flow
USE basic, ONLY: dt, chrono_ExBs, start_chrono, stop_chrono USE basic, ONLY: dt, time, chrono_ExBs, start_chrono, stop_chrono
USE grid, ONLY: local_nky, kyarray, inv_dkx USE grid, ONLY: local_nky, kyarray, inv_dkx, xarray,&
local_nkx, ikyarray, inv_ikyarray, deltakx, deltaky, deltax
USE model, ONLY: ExBrate USE model, ONLY: ExBrate
IMPLICIT NONE IMPLICIT NONE
! local var ! local var
INTEGER :: iky INTEGER :: iky, ix
REAL(xp):: dtExBshear
CALL start_chrono(chrono_ExBs) CALL start_chrono(chrono_ExBs)
! update the ExB shift, jumps and flags ! update the ExB shift, jumps and flags
shiftnow_ExB = .FALSE. shiftnow_ExB = .FALSE.
...@@ -125,6 +144,12 @@ CONTAINS ...@@ -125,6 +144,12 @@ CONTAINS
! in shiftnow_ExB and will use it in Shift_fields to avoid ! in shiftnow_ExB and will use it in Shift_fields to avoid
! zero-shiftings that may be majoritary. ! zero-shiftings that may be majoritary.
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0) shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
! Update the ExB nonlinear factor
dtExBshear = time - t0*inv_ikyarray(iky)*ANINT(ikyarray(iky)*time*inv_t0,xp)
DO ix = 1,local_nkx
ExB_NL_factor(ix,iky) = EXP(-imagu*xarray(ix)*ExBrate*ikyarray(iky)*dtExBshear)
inv_ExB_NL_factor(ix,iky) = 1._xp/ExB_NL_factor(ix,iky)
ENDDO
ENDDO ENDDO
CALL stop_chrono(chrono_ExBs) CALL stop_chrono(chrono_ExBs)
END SUBROUTINE Update_ExB_shear_flow END SUBROUTINE Update_ExB_shear_flow
......
src/fourier_mod.F90
+ 240
− 189
View file @ 1dccb6f8
This diff is collapsed.
src/grid_mod.F90
+ 18
− 4
View file @ 1dccb6f8
...@@ -25,7 +25,9 @@ MODULE grid ...@@ -25,7 +25,9 @@ MODULE grid
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: jarray, jarray_full INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: jarray, jarray_full
REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray ! ExB shear makes it ky dependant REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray ! ExB shear makes it ky dependant
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray_full REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray_full
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: xarray
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: ikyarray, inv_ikyarray !mode indices arrays
REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full
REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp
...@@ -72,7 +74,7 @@ MODULE grid ...@@ -72,7 +74,7 @@ MODULE grid
integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr_offset, local_nky_ptr_offset integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr_offset, local_nky_ptr_offset
! Grid spacing and limits ! Grid spacing and limits
REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz, inv_dkx REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz, inv_dkx
REAL(xp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max REAL(xp), PUBLIC, PROTECTED :: deltakx, deltaky, deltax, kx_max, ky_max, kx_min, ky_min!, kp_max
INTEGER , PUBLIC, PROTECTED :: local_pmin, local_pmax INTEGER , PUBLIC, PROTECTED :: local_pmin, local_pmax
INTEGER , PUBLIC, PROTECTED :: local_jmin, local_jmax INTEGER , PUBLIC, PROTECTED :: local_jmin, local_jmax
REAL(xp), PUBLIC, PROTECTED :: local_kymin, local_kymax REAL(xp), PUBLIC, PROTECTED :: local_kymin, local_kymax
...@@ -214,6 +216,8 @@ CONTAINS ...@@ -214,6 +216,8 @@ CONTAINS
local_nky_offset = local_nky_ptr_offset local_nky_offset = local_nky_ptr_offset
ALLOCATE(kyarray_full(Nky)) ALLOCATE(kyarray_full(Nky))
ALLOCATE(kyarray(local_nky)) ALLOCATE(kyarray(local_nky))
ALLOCATE(ikyarray(local_nky))
ALLOCATE(inv_ikyarray(local_nky))
ALLOCATE(AA_y(local_nky)) ALLOCATE(AA_y(local_nky))
!!---------------- RADIAL KX INDICES (not parallelized) !!---------------- RADIAL KX INDICES (not parallelized)
Nkx = Nx Nkx = Nx
...@@ -223,8 +227,9 @@ CONTAINS ...@@ -223,8 +227,9 @@ CONTAINS
local_nkx_ptr = ikxe - ikxs + 1 local_nkx_ptr = ikxe - ikxs + 1
local_nkx = ikxe - ikxs + 1 local_nkx = ikxe - ikxs + 1
local_nkx_offset = ikxs - 1 local_nkx_offset = ikxs - 1
ALLOCATE(kxarray(local_nky,local_Nkx))
ALLOCATE(kxarray_full(total_nkx)) ALLOCATE(kxarray_full(total_nkx))
ALLOCATE(kxarray(local_nky,local_Nkx))
ALLOCATE(xarray(total_nkx))
ALLOCATE(AA_x(local_nkx)) ALLOCATE(AA_x(local_nkx))
!!---------------- PARALLEL Z GRID (parallelized) !!---------------- PARALLEL Z GRID (parallelized)
total_nz = Nz total_nz = Nz
...@@ -432,10 +437,17 @@ CONTAINS ...@@ -432,10 +437,17 @@ CONTAINS
! indexation (|1 2 3||1 2 3|... local_nky|) ! indexation (|1 2 3||1 2 3|... local_nky|)
IF(Ny .EQ. 1) THEN IF(Ny .EQ. 1) THEN
kyarray(iky) = deltaky kyarray(iky) = deltaky
kyarray(iky) = iky-1
kyarray_full(iky) = deltaky kyarray_full(iky) = deltaky
SINGLE_KY = .TRUE. SINGLE_KY = .TRUE.
ELSE ELSE
kyarray(iky) = kyarray_full(iky+local_nky_offset) kyarray(iky) = kyarray_full(iky+local_nky_offset)
ikyarray(iky) = REAL((iky+local_nky_offset)-1,xp)
IF(ikyarray(iky) .GT. 0) THEN
inv_ikyarray(iky) = 1._xp/ikyarray(iky)
ELSE
inv_ikyarray(iky) = 0._xp
ENDIF
ENDIF ENDIF
! Finding kx=0 ! Finding kx=0
IF (kyarray(iky) .EQ. 0) THEN IF (kyarray(iky) .EQ. 0) THEN
...@@ -490,12 +502,14 @@ CONTAINS ...@@ -490,12 +502,14 @@ CONTAINS
Lx = Lx_adapted*Nexc Lx = Lx_adapted*Nexc
ENDIF ENDIF
deltakx = 2._xp*PI/Lx deltakx = 2._xp*PI/Lx
inv_dkx = 1._xp/deltakx inv_dkx = 1._xp/deltakx
deltax = Lx/REAL(Nx,xp) ! periodic donc pas Lx/(Nx-1)
IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3) IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
! Creating a grid ordered as dk*(0 1 2 3 -2 -1) ! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
DO ikx = 1,total_nkx DO ikx = 1,total_nkx
kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),xp) kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),xp)
IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx) IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
xarray(ikx) = REAL(ikx-1,xp)*deltax
END DO END DO
kx_max = MAXVAL(kxarray_full)!(total_nkx/2)*deltakx kx_max = MAXVAL(kxarray_full)!(total_nkx/2)*deltakx
kx_min = MINVAL(kxarray_full)!-kx_max+deltakx kx_min = MINVAL(kxarray_full)!-kx_max+deltakx
......
src/model_mod.F90
+ 5
− 0
View file @ 1dccb6f8
...@@ -31,6 +31,7 @@ MODULE model ...@@ -31,6 +31,7 @@ MODULE model
! Auxiliary variable ! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift
LOGICAL, PUBLIC, PROTECTED :: ExB = .false. ! presence of ExB background shearing rate
! Removes Landau damping in temperature and higher equation (Ivanov 2022) ! Removes Landau damping in temperature and higher equation (Ivanov 2022)
LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false. LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false.
! Flag to force the reality condition symmetry for the kx at ky=0 ! Flag to force the reality condition symmetry for the kx at ky=0
...@@ -78,6 +79,10 @@ CONTAINS ...@@ -78,6 +79,10 @@ CONTAINS
EM = .FALSE. EM = .FALSE.
ENDIF ENDIF
IF(ExBrate .GT. 0) THEN
ExB = .TRUE.
ENDIF
END SUBROUTINE model_readinputs END SUBROUTINE model_readinputs
SUBROUTINE model_outputinputs(fid) SUBROUTINE model_outputinputs(fid)
......
src/nonlinear_mod.F90
+ 15
− 14
View file @ 1dccb6f8
MODULE nonlinear MODULE nonlinear
USE array, ONLY : dnjs, Sapj, kernel USE array, ONLY : dnjs, Sapj, kernel
USE fourier, ONLY : bracket_sum_r, bracket_sum_c, planf, planb, poisson_bracket_and_sum USE fourier, ONLY : bracket_sum_r, bracket_sum_c, planf, planb, poisson_bracket_and_sum,&
apply_inv_ExB_NL_factor
USE fields, ONLY : phi, psi, moments USE fields, ONLY : phi, psi, moments
USE grid, ONLY: local_na, & USE grid, ONLY : local_na, &
local_np,ngp,parray,pmax,& local_np,ngp,parray,pmax,&
local_nj,ngj,jarray,jmax, local_nj_offset, dmax,& local_nj,ngj,jarray,jmax, local_nj_offset, dmax,&
kyarray, AA_y, local_nky_ptr, local_nky_ptr_offset,inv_Ny,& kyarray, AA_y, local_nky_ptr, local_nky_ptr_offset,inv_Ny,&
local_nkx_ptr,kxarray, AA_x, inv_Nx,& local_nkx_ptr,kxarray, AA_x, inv_Nx,&
local_nz,ngz,zarray,nzgrid, deltakx, iky0, contains_kx0, contains_ky0 local_nz,ngz,zarray,nzgrid, deltakx, iky0, contains_kx0, contains_ky0
USE model, ONLY : LINEARITY, EM, ikxZF, ZFamp USE model, ONLY : LINEARITY, EM, ikxZF, ZFamp, ExB
USE closure, ONLY : evolve_mom, nmaxarray USE closure, ONLY : evolve_mom, nmaxarray
USE prec_const, ONLY : xp USE prec_const, ONLY : xp
USE species, ONLY : sqrt_tau_o_sigma USE species, ONLY : sqrt_tau_o_sigma
USE time_integration, ONLY : updatetlevel USE time_integration, ONLY : updatetlevel
USE ExB_shear_flow, ONLY : ExB_NL_factor, inv_ExB_NL_factor
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
IMPLICIT NONE IMPLICIT NONE
...@@ -20,8 +22,6 @@ MODULE nonlinear ...@@ -20,8 +22,6 @@ MODULE nonlinear
INCLUDE 'fftw3-mpi.f03' INCLUDE 'fftw3-mpi.f03'
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: Fx_cmpx, Gy_cmpx
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: Fy_cmpx, Gx_cmpx, F_conv_G
INTEGER :: in, is, p_int, j_int, n_int INTEGER :: in, is, p_int, j_int, n_int
INTEGER :: smax INTEGER :: smax
REAL(xp):: sqrt_p, sqrt_pp1 REAL(xp):: sqrt_p, sqrt_pp1
...@@ -33,13 +33,6 @@ SUBROUTINE nonlinear_init ...@@ -33,13 +33,6 @@ SUBROUTINE nonlinear_init
IMPLICIT NONE IMPLICIT NONE
ALLOCATE( F_cmpx(local_nky_ptr,local_nkx_ptr)) ALLOCATE( F_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE( G_cmpx(local_nky_ptr,local_nkx_ptr)) ALLOCATE( G_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE(Fx_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE(Gy_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE(Fy_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE(Gx_cmpx(local_nky_ptr,local_nkx_ptr))
ALLOCATE(F_conv_G(local_nky_ptr,local_nkx_ptr))
END SUBROUTINE nonlinear_init END SUBROUTINE nonlinear_init
SUBROUTINE compute_Sapj SUBROUTINE compute_Sapj
...@@ -105,7 +98,9 @@ SUBROUTINE compute_nonlinear ...@@ -105,7 +98,9 @@ SUBROUTINE compute_nonlinear
dnjs(in,ij,is) * moments(ia,ipi,isi,:,:,izi,updatetlevel) dnjs(in,ij,is) * moments(ia,ipi,isi,:,:,izi,updatetlevel)
ENDDO s1 ENDDO s1
! this function adds its result to bracket_sum_r ! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r) CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,&
ExB, ExB_NL_factor, bracket_sum_r)
!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi} !-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
IF(EM) THEN IF(EM) THEN
! First convolution terms ! First convolution terms
...@@ -119,9 +114,15 @@ SUBROUTINE compute_nonlinear ...@@ -119,9 +114,15 @@ SUBROUTINE compute_nonlinear
+sqrt_p *moments(ia,ipi-1,isi,:,:,izi,updatetlevel)) +sqrt_p *moments(ia,ipi-1,isi,:,:,izi,updatetlevel))
ENDDO s2 ENDDO s2
! this function adds its result to bracket_sum_r ! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r) CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,&
ExB, ExB_NL_factor,bracket_sum_r)
ENDIF ENDIF
ENDDO n ENDDO n
! Apply the ExB shearing rate factor before going back to k-space
IF (ExB) THEN
CALL apply_inv_ExB_NL_factor(bracket_sum_r,inv_ExB_NL_factor)
ENDIF
! Put the real nonlinear product back into k-space ! Put the real nonlinear product back into k-space
#ifdef SINGLE_PRECISION #ifdef SINGLE_PRECISION
call fftwf_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) call fftwf_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment