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

save current changes (ExB still bug)

parent 8270b2c6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ExB_shear_flow_mod.F90
+ 34
− 24
View file @ e37f439b
......@@ -9,6 +9,7 @@ MODULE ExB_shear_flow
REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E
REAL(xp), PUBLIC, PROTECTED :: t0, inv_t0 = 0._xp ! charact. shear time
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky)
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB_full ! full ky version
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: dkx_ExB ! correction to obtain the exact kx mode
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
......@@ -16,19 +17,20 @@ MODULE ExB_shear_flow
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: inv_ExB_NL_factor
LOGICAL, PUBLIC, PROTECTED :: ExB = .false. ! presence of ExB background shearing rate
! Routines
PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
PUBLIC :: Setup_ExB_shear_flow, Array_shift_ExB_shear_flow, Update_ExB_shear_flow
CONTAINS
! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow(ExBrate)
USE grid, ONLY : Nx, local_nky, local_nx, Ny, deltakx, deltaky
USE grid, ONLY: Nx, local_nky, total_nky, local_nx, Ny, deltakx, deltaky
USE geometry, ONLY: Cyq0_x0
IMPLICIT NONE
REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shearing rate and aux var
gamma_E = ExBrate
IF(gamma_E .GT. 0._xp) THEN
gamma_E = -ExBrate*Cyq0_x0
IF(abs(gamma_E) .GT. EPSILON(gamma_E)) THEN
ExB = .TRUE.
t0 = deltakx/deltaky/gamma_E
inv_t0 = 1._xp/t0
......@@ -41,6 +43,8 @@ CONTAINS
! Setup the ExB shift array
ALLOCATE(sky_ExB(local_nky))
sky_ExB = 0._xp
ALLOCATE(sky_ExB_full(total_nky))
sky_ExB_full = 0._xp
! Setup the kx correction array
ALLOCATE(dkx_ExB(local_nky))
......@@ -63,16 +67,14 @@ CONTAINS
END SUBROUTINE Setup_ExB_shear_flow
! Update according to the current ExB shear value
! -the grids
! -shift the grids
! -the spatial operators
! -the fields by imposing a shift on kx
SUBROUTINE Apply_ExB_shear_flow
USE basic, ONLY: time
USE grid, ONLY: local_nky, kxarray, update_grids, &
total_nkx, deltakx, kx_min, kx_max
SUBROUTINE Array_shift_ExB_shear_flow
USE grid, ONLY: local_nky, update_grids, &
total_nkx, deltakx, kx_min, kx_max, kxarray0
USE prec_const, ONLY: PI
USE fields, ONLY: moments, phi, psi
USE geometry, ONLY: gxx,gxy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
IMPLICIT NONE
! local var
......@@ -105,10 +107,17 @@ CONTAINS
ikx = i_ - total_nkx/2 + 1
ENDIF
ikx_s = ikx + jump_ExB(iky)
! adjust the shift according
IF (ikx_s .LE. 0) &
ikx_s = ikx_s + total_nkx
IF (ikx_s .GT. total_nkx) &
ikx_s = ikx_s - total_nkx
! We test if the shifted modes are still in contained in our resolution
IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
(((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
! IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
! (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
! ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. kx_max) .AND. &
(kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. kx_min)) THEN
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
phi(iky,ikx,:) = phi(iky,ikx_s,:)
psi(iky,ikx,:) = psi(iky,ikx_s,:)
......@@ -123,20 +132,20 @@ CONTAINS
ENDIF
ENDDO
ENDIF
END SUBROUTINE Apply_ExB_shear_flow
END SUBROUTINE Array_shift_ExB_shear_flow
! update the ExB shear value for the next time step
SUBROUTINE Update_ExB_shear_flow
USE basic, ONLY: dt, time, chrono_ExBs, start_chrono, stop_chrono
USE grid, ONLY: local_nky, local_nky_offset, kyarray, kyarray_full, inv_dkx, xarray, Nx, Ny, &
local_nx, local_nx_offset, deltax, kxarray, &
ikyarray, inv_ikyarray, deltakx, deltaky, update_grids
USE basic, ONLY: dt, time
USE grid, ONLY: local_nky, local_nky_offset, kyarray, inv_dkx, xarray, Nx, Ny, &
local_nx, local_nx_offset, kyarray_full,&
ikyarray, inv_ikyarray, deltaky, update_grids
USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
IMPLICIT NONE
! local var
INTEGER :: iky, ix
REAL(xp):: dt_ExB, ky, kx, J_dp, inv_J, x, dtshift
REAL(xp):: dt_ExB, J_dp, inv_J, x
! do nothing if no ExB
IF(ExB) THEN
! reset the ExB shift, jumps and flags
......@@ -149,7 +158,8 @@ CONTAINS
! zero-shiftings that may be majoritary.
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
ENDDO
! Update the full skyExB array too
sky_ExB_full = sky_ExB_full - kyarray_full*gamma_E*dt
! Update the ExB nonlinear factor...
DO iky = 1,local_Nky ! WARNING: Local indices ky loop
! for readability
......@@ -161,7 +171,8 @@ CONTAINS
DO ix = 1,Nx
x = xarray(ix)
! assemble the ExB nonlin factor
ExB_NL_factor(ix,iky) = EXP(imagu*x*dkx_ExB(iky))
!ExB_NL_factor(ix,iky) = EXP(imagu*x*dkx_ExB(iky))
ExB_NL_factor(ix,iky) = EXP(imagu*x*sky_ExB(iky)) !GENE does that???
ENDDO
ENDDO
! ... and the inverse
......@@ -174,15 +185,14 @@ CONTAINS
DO ix = 1,local_nx
x = xarray(ix+local_nx_offset)
! assemble the inverse ExB nonlin factor
inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*gamma_E*J_dp*deltaky*dt_ExB)
! inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*gamma_E*J_dp*deltaky*dt_ExB)
inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*sky_ExB_full(iky)) !GENE does that???
ENDDO
ENDDO
ENDIF
! We update the operators and grids
! update the grids
CALL update_grids(dkx_ExB,gxx,gxy,gyy,inv_hatB2)
!write(42,*), kxarray(5,3)
! update the EM op., the kernels and the curvature op.
CALL evaluate_kernels
CALL evaluate_EM_op
......
src/auxval.F90
+ 2
− 2
View file @ e37f439b
......@@ -2,9 +2,9 @@ subroutine auxval
! Set auxiliary values, at beginning of simulation
USE basic, ONLY: str, speak
USE grid, ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset, total_nj,&
USE grid, ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset,&
local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset, dmax,&
local_nz, local_nz_offset, total_nz, init_grids_data, set_grids
local_nz, local_nz_offset, init_grids_data, set_grids
!USE array
USE model, ONLY: Na, EM, LINEARITY, N_HD, ExBrate
USE fourier, ONLY: init_grid_distr_and_plans
......
src/control.F90
+ 21
− 7
View file @ e37f439b
SUBROUTINE control
! Control the run
use basic, ONLY: str,daytim,speak,basic_data,&
use basic, ONLY: str,daytim,speak,basic_data,&
nlend,step,increase_step,increase_time,increase_cstep,&
chrono_runt,chrono_step, chrono_diag, start_chrono, stop_chrono
use prec_const, ONLY: xp, stdout
USE parallel, ONLY: ppinit
USE initial, ONLY: initialize
USE mpi, ONLY: MPI_COMM_WORLD
USE diagnostics, ONLY: diagnose
chrono_runt,chrono_step, chrono_diag, chrono_ExBs,&
start_chrono, stop_chrono
use prec_const, ONLY: xp, stdout
USE parallel, ONLY: ppinit
USE initial, ONLY: initialize
USE mpi, ONLY: MPI_COMM_WORLD
USE diagnostics, ONLY: diagnose
USE ExB_shear_flow, ONLY: Array_shift_ExB_shear_flow, Update_ExB_shear_flow
IMPLICIT NONE
REAL(xp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
......@@ -70,6 +72,18 @@ SUBROUTINE control
CALL tesend
IF( nlend ) EXIT ! exit do loop
! Update the ExB shear flow for the next step
! This call includes :
! - the ExB shear value (s(ky)) update for the next time step
! - the kx grid update
! - the ExB NL correction factor update (exp(+/- ixkySdts))
! - (optional) the kernel, poisson op. and ampere op update
CALL start_chrono(chrono_ExBs)
CALL Update_ExB_shear_flow
! Shift the arrays if the shear value sky is too high
CALL Array_shift_ExB_shear_flow
CALL stop_chrono(chrono_ExBs)
! Increment steps and csteps (private in basic module)
CALL increase_step
CALL increase_cstep
......
src/diagnostics_mod.F90
+ 2
− 4
View file @ e37f439b
......@@ -120,7 +120,7 @@ CONTAINS
SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
USE basic, ONLY: lu_in, chrono_runt, cstep, dt, time, tmax, display_h_min_s
USE basic, ONLY: lu_in, chrono_runt, cstep, time, display_h_min_s
USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id
IMPLICIT NONE
......@@ -159,7 +159,7 @@ CONTAINS
USE basic, ONLY: speak,chrono_runt,&
cstep,iframe0d,iframe3d,iframe5d,crashed
USE grid, ONLY: &
parray_full,pmax,jarray_full,jmax, kparray, &
parray_full,jarray_full, kparray, &
kyarray_full,kxarray_full,zarray_full, ngz, total_nz, local_nz, ieven,&
local_Nky, total_nky, local_nkx, total_nkx
USE geometry, ONLY: gxx, gxy, gxz, gyy, gyz, gzz, &
......@@ -390,7 +390,6 @@ CONTAINS
USE futils, ONLY: append, attach, getatt, creatd
USE prec_const
USE processing
USE model, ONLY: Na
IMPLICIT NONE
! Time measurement data
CALL append(fidres, "/profiler/Tc_rhs", REAL(chrono_mrhs%ttot,dp),ionode=0)
......@@ -436,7 +435,6 @@ CONTAINS
USE CLA, ONLY: Sf
#endif
IMPLICIT NONE
CHARACTER(50) :: dset_name
iframe2d=iframe2d+1
CALL append(fidres,"/data/var2d/time", REAL(time,dp), ionode=0)
CALL append(fidres,"/data/var2d/cstep",REAL(cstep,dp),ionode=0)
......
src/fourier_mod.F90
+ 59
− 33
View file @ e37f439b
......@@ -44,11 +44,10 @@ MODULE fourier
!******************************************************************************!
!------------- Initialize the grid distribution and plans -------------
! If we use the 2D fftw routines, the fftw library decides the best data distribution
SUBROUTINE init_grid_distr_and_plans(FFT2D,Nx,Ny,communicator,&
SUBROUTINE init_grid_distr_and_plans(Nx,Ny,communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
USE basic, ONLY: speak
IMPLICIT NONE
LOGICAL, INTENT(IN) :: FFT2D
INTEGER, INTENT(IN) :: Nx,Ny, communicator
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nx_ptr,local_nx_ptr_offset
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nky_ptr,local_nky_ptr_offset
......@@ -57,7 +56,7 @@ MODULE fourier
NY_halved = NY_/2 + 1
inv_Ny_ = 1._xp/REAL(2*NY_halved,xp)
! Call FFTW 2D mpi routines to distribute the data and init 2D MPI FFTW plans
CALL fft2D_distr_and_plans(Nx,Ny,communicator,&
CALL fft2D_distr_and_plans(communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
local_nx_ = local_nx_ptr ! store number of local x in the module
local_nky_ = local_nky_ptr ! store number of local ky in the module
......@@ -67,10 +66,10 @@ MODULE fourier
END SUBROUTINE init_grid_distr_and_plans
!------------- 2D fft initialization and mpi distribution
SUBROUTINE fft2D_distr_and_plans(Nx,Ny,communicator,&
SUBROUTINE fft2D_distr_and_plans(communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
IMPLICIT NONE
INTEGER, INTENT(IN) :: Nx,Ny, communicator
INTEGER, INTENT(IN) :: communicator
! Distribution in the real space along x
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nx_ptr, local_nx_ptr_offset
! Distribution in the fourier space along ky
......@@ -138,18 +137,14 @@ MODULE fourier
!------------- 1D initialization with balanced data distribution
SUBROUTINE fft1D_plans
USE utility, ONLY: decomp1D
USE parallel, ONLY: num_procs_ky, rank_ky
IMPLICIT NONE
! local var
integer(C_INTPTR_T) :: rank ! rank of each 1D fourier transforms
integer(C_INTPTR_T) :: n ! size of the data to fft
integer(C_INTPTR_T) :: howmany ! howmany 1D fourier transforms
COMPLEX, DIMENSION(:,:), ALLOCATABLE:: in, out
integer(C_INTPTR_T) :: inembed, onembed
integer(C_INTPTR_T) :: istride, ostride
integer(C_INTPTR_T) :: idist, odist
integer(C_INTPTR_T) :: sign
integer(C_INTPTR_T) :: flags
!! Plan of the 4 1D many transforms required
!----------- 1: FFTx and inv through local ky data
......@@ -273,12 +268,11 @@ END SUBROUTINE fft1D_plans
! module variable (convolution theorem)
SUBROUTINE poisson_bracket_and_sum( ky_array, kx_array, inv_Ny, inv_Nx, AA_y, AA_x,&
local_nky, total_nkx, F_, G_,&
ExB, ExB_NL_factor, dkx_ExB, sum_real_)
USE parallel, ONLY: my_id, num_procs_ky, comm_ky, rank_ky
ExB, ExB_NL_factor, sum_real_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: local_nky,total_nkx
REAL(xp), INTENT(IN) :: inv_Nx, inv_Ny
REAL(xp), DIMENSION(local_nky), INTENT(IN) :: ky_array, AA_y, dkx_ExB
REAL(xp), DIMENSION(local_nky), INTENT(IN) :: ky_array, AA_y
REAL(xp), DIMENSION(total_nkx), INTENT(IN) :: AA_x
REAL(xp), DIMENSION(local_nky,total_nkx), INTENT(IN) :: kx_array
COMPLEX(c_xp_c), DIMENSION(local_nky,total_nkx), &
......@@ -298,10 +292,10 @@ END SUBROUTINE fft1D_plans
DO ikx = 1,total_nkx
ky = ky_array(iky)
kxs = kx_array(iky,ikx)
ikxF(ikx,iky) = imagu*kxs*F_(iky,ikx)*AA_y(iky)*AA_x(ikx)
ikyG(ikx,iky) = imagu*ky *G_(iky,ikx)*AA_y(iky)*AA_x(ikx)
ikyF(ikx,iky) = imagu*ky *F_(iky,ikx)*AA_y(iky)*AA_x(ikx)
ikxG(ikx,iky) = imagu*kxs*G_(iky,ikx)*AA_y(iky)*AA_x(ikx)
ikxF(ikx,iky) = imagu*kxs*F_(iky,ikx)
ikyG(ikx,iky) = imagu*ky *G_(iky,ikx)
ikyF(ikx,iky) = imagu*ky *F_(iky,ikx)
ikxG(ikx,iky) = imagu*kxs*G_(iky,ikx)
ENDDO
ENDDO
IF(ExB) THEN
......@@ -311,6 +305,13 @@ END SUBROUTINE fft1D_plans
CALL apply_ExB_NL_factor(ikyF,ExB_NL_factor)
CALL apply_ExB_NL_factor(ikxG,ExB_NL_factor)
ENDIF
! Anti Aliasing
DO iky = 1,local_nky
ikxF(:,iky) = ikxF(:,iky)*AA_y(iky)*AA_x(:)
ikyG(:,iky) = ikyG(:,iky)*AA_y(iky)*AA_x(:)
ikyF(:,iky) = ikyF(:,iky)*AA_y(iky)*AA_x(:)
ikxG(:,iky) = ikxG(:,iky)*AA_y(iky)*AA_x(:)
ENDDO
!-------------------- First term df/dx x dg/dy --------------------
#ifdef SINGLE_PRECISION
call fftwf_mpi_execute_dft_c2r(planb, ikxF, real_data_f)
......@@ -342,38 +343,63 @@ END SUBROUTINE fft1D_plans
COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(INOUT) :: fkxky
COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: ExB_NL_factor
! local variables
COMPLEX(xp), DIMENSION(NX_,local_nky_) :: fxky
CALL iFFT_kxky_to_xky(fkxky,fxky)
fxky = fxky*ExB_NL_factor*inv_Nx_
CALL FFT_xky_to_kxky(fxky,fkxky)
COMPLEX(xp), DIMENSION(NX_,local_nky_) :: tmp_kxky, tmp_xky
INTEGER :: ix,ikx,iky
! Fill the buffer
DO iky = 1,local_nky_
DO ikx = 1,NX_
tmp_kxky(ikx,iky) = fkxky(ikx,iky)
ENDDO
ENDDO
! go to x,ky
CALL iFFT_kxky_to_xky(tmp_kxky,tmp_xky)
! multiply result by NL factor
DO iky = 1,local_nky_
DO ix = 1,NX_
tmp_xky(ix,iky) = inv_Nx_*tmp_xky(ix,iky)*ExB_NL_factor(ix,iky)
ENDDO
ENDDO
! back to kx,ky
CALL FFT_xky_to_kxky(tmp_xky,tmp_kxky)
! fill output array
DO iky = 1,local_nky_
DO ikx = 1,NX_
fkxky(ikx,iky) = tmp_kxky(ikx,iky)
ENDDO
ENDDO
END SUBROUTINE apply_ExB_NL_factor
! Apply the exp(i*x*ky*S*J*dt) factor to the real Poisson Bracket sum [f,g]
SUBROUTINE apply_inv_ExB_NL_factor(fyx,inv_ExB_NL_factor)
IMPLICIT NONE
real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
COMPLEX(xp), DIMENSION(NY_halved,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
! local variables
REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: buffer_1, buffer_2
COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: fkyx
INTEGER :: ix, iy
! print*, SIZE(fyx), SIZE(buffer_1)
REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: tmp_yx_1, tmp_yx_2
COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: tmp_kyx
INTEGER :: ix, iy, iky
! Fill buffer
DO ix = 1,local_nx_
DO iy = 1,2*NY_halved
buffer_1(iy,ix) = fyx(iy,ix)
tmp_yx_1(iy,ix) = fyx(iy,ix)
ENDDO
ENDDO
! print*, SUM(fyx), SUM(buffer_1)
CALL FFT_yx_to_kyx(buffer_1,fkyx)
fkyx = fkyx*inv_ExB_NL_factor
CALL iFFT_kyx_to_yx(fkyx,buffer_2)
! Fourier real to complex on the second buffer (first buffer is now unusable)
CALL FFT_yx_to_kyx(tmp_yx_1,tmp_kyx)
! Treat the result with the ExB NL factor
DO iky = 1,NY_halved
DO ix = 1,local_nx_
tmp_kyx(iky,ix) = tmp_kyx(iky,ix)*inv_ExB_NL_factor(iky,ix)
ENDDO
ENDDO
! Back to Fourier space in the third buffer
CALL iFFT_kyx_to_yx(tmp_kyx,tmp_yx_2)
! Fill the result array with normalization of iFFT
DO ix = 1,local_nx_
DO iy = 1,2*NY_halved
fyx(iy,ix) = buffer_2(iy,ix)*inv_Ny_
fyx(iy,ix) = tmp_yx_2(iy,ix)*inv_Ny_
ENDDO
ENDDO
! print*, SUM(fyx)
! stop
END SUBROUTINE apply_inv_ExB_NL_factor
!******************************************************************************!
......
src/grid_mod.F90
+ 7
− 13
View file @ e37f439b
......@@ -112,11 +112,6 @@ MODULE grid
! 1D Antialiasing arrays (2/3 rule)
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_x
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_y
! 2D fft routines or 1D methods (for ExBshear simulations, 1D is required)
! The 2D fft is using fftw mpi optimization
! The 1D is not using mpi and does a data transfer with redundant computations
! (each process computes the convolution)
LOGICAL, PUBLIC, PROTECTED :: FFT2D = .TRUE.
! Public Functions
PUBLIC :: init_grids_data, set_grids, update_grids
......@@ -134,7 +129,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER :: lun = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmax, jmax, Nx, Lx, Ny, Ly, Nz, SG, Nexc, FFT2D
NAMELIST /GRID/ pmax, jmax, Nx, Lx, Ny, Ly, Nz, SG, Nexc
READ(lun,grid)
IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use
SG = .FALSE.
......@@ -202,7 +197,7 @@ CONTAINS
!! Parallel distribution of kx ky grid
IF (LINEARITY .NE. 'linear') THEN ! we let FFTW distribute if we use it
IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
CALL init_grid_distr_and_plans(FFT2D,Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
CALL init_grid_distr_and_plans(Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
ELSE ! otherwise we distribute equally
IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
! balanced distribution among the processes
......@@ -453,7 +448,6 @@ CONTAINS
! indexation (|1 2 3||1 2 3|... local_nky|)
IF(Ny .EQ. 1) THEN
kyarray(iky) = deltaky
kyarray(iky) = iky-1
kyarray_full(iky) = deltaky
SINGLE_KY = .TRUE.
ELSE
......@@ -543,7 +537,7 @@ CONTAINS
ERROR STOP "Gyacomo is safer with an even Kx number"
ENDIF
! Orszag 2/3 filter
two_third_kxmax = 2._xp/3._xp*kx_max;
two_third_kxmax = 2._xp/3._xp*(kx_max-deltakx);
! Antialiasing filter
DO iky = 1,local_nky
DO ikx = 1,local_nkx
......@@ -673,11 +667,11 @@ CONTAINS
REAL(xp), DIMENSION(local_nky),INTENT(IN) :: dkx_ExB ! ExB correction dkx = gamma_E ky dtshift
REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
INTEGER :: eo,iz,iky,ikx
REAL(xp) :: kx, ky, skp2
REAL(xp) :: kx, ky
! Update the kx grid
DO iky = 1,local_nky
DO ikx = 1,total_Nkx
kxarray(iky,ikx) = kxarray0(ikx) + dkx_ExB(iky)
DO ikx = 1,total_Nkx
DO iky = 1,local_nky
kxarray(iky,ikx) = kxarray0(ikx) - dkx_ExB(iky)
ENDDO
ENDDO
! Update the kperp grid
......
src/initial_mod.F90
+ 3
− 13
View file @ e37f439b
......@@ -46,7 +46,6 @@ CONTAINS
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
USE restarts, ONLY: load_moments, job2load
USE processing, ONLY: compute_fluid_moments, compute_radial_heatflux, compute_radial_transport
USE model, ONLY: LINEARITY
USE nonlinear, ONLY: compute_Sapj, nonlinear_init
IMPLICIT NONE
CALL set_updatetlevel(1)
......@@ -266,19 +265,10 @@ CONTAINS
!!!!!!! Initialize a single mode in the gyrocenter density
!******************************************************************************!
SUBROUTINE init_single_mode
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
ngp, ngj, ngz, iky0, kxarray, kyarray, parray, jarray, &
contains_ky0, deltakx, deltaky
USE fields, ONLY: moments
USE prec_const, ONLY: xp
USE model, ONLY: LINEARITY
USE parallel, ONLY: my_id
IMPLICIT NONE
REAL(xp) :: amplitude
INTEGER :: ia,ip,ij,ikx,iky,iz, ikxi, ikyi
INTEGER, DIMENSION(12) :: iseedarr
moments = 0._xp
IF (my_id .EQ. 0) THEN
moments(:,:,:,2,2,:,:) = init_amp
......@@ -470,7 +460,7 @@ CONTAINS
USE fields, ONLY: moments
USE prec_const, ONLY: xp, pi
USE model, ONLY: LINEARITY
USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
USE geometry, ONLY: Jacobian, iInt_Jacobian
IMPLICIT NONE
REAL(xp) :: kx, ky, sigma_z, z
INTEGER :: ia,iky,ikx,iz,ip,ij
......@@ -548,8 +538,8 @@ CONTAINS
USE basic, ONLY: maindir
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
local_nkx_offset, local_nky_offset, kxarray, kyarray, &
ngp, ngj, ngz, iky0, ikx0, parray, jarray,&
deltakx, deltaky, contains_ky0, contains_kx0,&
ngp, ngj, ngz, parray, jarray,&
deltakx, deltaky,&
AA_x, AA_y
USE fields, ONLY: moments
USE prec_const, ONLY: xp, imagu
......
src/nonlinear_mod.F90
+ 2
− 2
View file @ e37f439b
......@@ -100,7 +100,7 @@ SUBROUTINE compute_nonlinear
! 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,total_nkx,F_cmpx,G_cmpx,&
ExB, ExB_NL_factor, dkx_ExB, bracket_sum_r)
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}
IF(EM) THEN
! First convolution terms
......@@ -116,7 +116,7 @@ SUBROUTINE compute_nonlinear
! 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,total_nkx,F_cmpx,G_cmpx,&
ExB, ExB_NL_factor, dkx_ExB, bracket_sum_r)
ExB, ExB_NL_factor, bracket_sum_r)
ENDIF
ENDDO n
! Apply the ExB shearing rate factor before going back to k-space
......
src/numerics_mod.F90
+ 1
− 0
View file @ e37f439b
......@@ -63,6 +63,7 @@ SUBROUTINE build_dv4Hp_table
!we scale it w.r.t. to the max degree since
!D_4^{v}\sim (\Delta v/2)^4 and \Delta v \sim 2pi/kvpar = pi/\sqrt{2P}
! dv4_Hp_coeff = dv4_Hp_coeff*(1._xp/2._xp/SQRT(REAL(pmax,xp)))**4
IF(pmax .GT. 0) &
dv4_Hp_coeff = dv4_Hp_coeff*(PI/2._xp/SQRT(2._xp*REAL(pmax,xp)))**4
END SUBROUTINE build_dv4Hp_table
!******************************************************************************!
......
src/stepon.F90
+ 3
− 19
View file @ e37f439b
......@@ -2,14 +2,13 @@ SUBROUTINE stepon
! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
USE advance_field_routine, ONLY: advance_time_level, advance_moments
USE basic, ONLY: nlend, chrono_advf, chrono_pois,&
chrono_chck, chrono_clos, chrono_ghst, chrono_ExBs,&
chrono_chck, chrono_clos, chrono_ghst,&
start_chrono, stop_chrono
USE closure, ONLY: apply_closure_model
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: xp
USE ExB_shear_flow, ONLY: Apply_ExB_shear_flow, Update_ExB_shear_flow
#ifdef TEST_SVD
USE CLA, ONLY: test_svd,filter_sv_moments_ky_pj
#endif
......@@ -18,17 +17,7 @@ SUBROUTINE stepon
INTEGER :: num_step, ierr
LOGICAL :: mlend
! Update the ExB backg. shear flow before the step
! This call includes :
! - the kx grid update
! - the kernel, poisson op. and ampere op update
! - kx-shift of the fields at required ky values
! - the ExB shear value (s(ky)) update for the next time step
CALL start_chrono(chrono_ExBs)
CALL Apply_ExB_shear_flow
CALL stop_chrono(chrono_ExBs)
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
SUBSTEPS:DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
......@@ -80,12 +69,7 @@ SUBROUTINE stepon
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
!----- AFTER: All fields are updated for step = n+1
END DO
! Update the ExB shear flow for the next step
CALL start_chrono(chrono_ExBs)
CALL Update_ExB_shear_flow
CALL stop_chrono(chrono_ExBs)
END DO SUBSTEPS
CONTAINS
!!!! Basic structure to simplify stepon
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment