Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (half_ky_plane) Browse

gyacomo/src/inital.F90
96dc3536888dhalf_ky_plane

Last Change
Raw File

inital.F90
View Options

!******************************************************************************!
!!!!!! initialize the moments and load/build coeff tables
!******************************************************************************!
SUBROUTINE inital
USE basic
USE model, ONLY : CO, NON_LIN
USE initial_par
USE prec_const
USE coeff
USE time_integration
USE array, ONLY : Sepj,Sipj
USE collision
USE closure
USE ghosts
USE restarts
implicit none
CALL set_updatetlevel(1)
IF (my_id .EQ. 0) WRITE(*,*) 'Evaluate kernels'
CALL evaluate_kernels
!!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
IF ( RESTART ) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
CALL load_moments ! get N_0
IF (my_id .EQ. 0) WRITE(*,*) 'Init phi with Poisson'
CALL poisson ! compute phi_0=phi(N_0)
ELSE
IF (INIT_NOISY_PHI) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi ! init noisy phi_0, N_0 = 0
ELSEIF (INIT_ZF .GT. 0) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
CALL init_phi ! init ZF
ELSE
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments and ghosts'
CALL init_moments ! init noisy N_0
IF (my_id .EQ. 0) WRITE(*,*) 'Init phi with Poisson'
CALL poisson ! get phi_0 = phi(N_0)
ENDIF
ENDIF
IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
CALL apply_closure_model
IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
CALL update_ghosts
!!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
IF ( NON_LIN ) THEN;
IF (my_id .EQ. 0) WRITE(*,*) 'Building Dnjs table'
CALL build_dnjs_table
IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
ENDIF
!!!!!! Load the COSOlver collision operator coefficients !!!!!!
IF (ABS(CO) .GT. 1) THEN
CALL load_COSOlver_mat
! Compute collision
CALL compute_TColl ! compute C_0 = C(N_0)
ENDIF
END SUBROUTINE inital
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize the moments randomly
!******************************************************************************!
SUBROUTINE init_moments
USE basic
USE grid
USE fields
USE prec_const
USE utility, ONLY: checkfield
USE initial_par
USE model, ONLY : NON_LIN
IMPLICIT NONE
REAL(dp) :: noise
REAL(dp) :: kr, kz, sigma, gain, kz_shift
INTEGER, DIMENSION(12) :: iseedarr
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr+my_id)
!**** Broad noise initialization *******************************************
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
moments_e( ip,ij, ikr,ikz, :) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
IF ( contains_kr0 ) THEN
DO ikz=2,Nkz/2 !symmetry at kr = 0
moments_e( ip,ij,ikr_0,ikz, :) = moments_e( ip,ij,ikr_0,Nkz+2-ikz, :)
END DO
ENDIF
END DO
END DO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
moments_i( ip,ij, ikr,ikz, :) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
IF ( contains_kr0 ) THEN
DO ikz=2,Nkz/2 !symmetry at kr = 0
moments_i( ip,ij,ikr_0,ikz, :) = moments_i( ip,ij,ikr_0,Nkz+2-ikz, :)
END DO
ENDIF
END DO
END DO
! Putting to zero modes that are not in the 2/3 Orszag rule
IF (NON_LIN) THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
moments_e( ip,ij,ikr,ikz, :) = moments_e( ip,ij,ikr,ikz, :)*AA_r(ikr)*AA_z(ikz)
ENDDO
ENDDO
ENDDO
ENDDO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
moments_i( ip,ij,ikr,ikz, :) = moments_i( ip,ij,ikr,ikz, :)*AA_r(ikr)*AA_z(ikz)
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE init_moments
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize a noisy ES potential and cancel the moments
!******************************************************************************!
SUBROUTINE init_phi
USE basic
USE grid
USE fields
USE prec_const
USE initial_par
IMPLICIT NONE
REAL(dp) :: noise
REAL(dp) :: kr, kz, sigma, gain, kz_shift
INTEGER, DIMENSION(12) :: iseedarr
IF (INIT_NOISY_PHI) THEN
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr+my_id)
!**** noise initialization *******************************************
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
phi(ikr,ikz) = (init_background + init_noiselvl*(noise-0.5_dp))*AA_r(ikr)*AA_z(ikz)
END DO
END DO
!symmetry at kr = 0 to keep real inverse transform
IF ( contains_kr0 ) THEN
DO ikz=2,Nkz/2
phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz)
END DO
phi(ikr_0,Nz/2) = REAL(phi(ikr_0,Nz/2)) !origin must be real
ENDIF
!**** Cancel previous moments initialization
moments_e = 0._dp; moments_i = 0._dp
ELSEIF(INIT_ZF .GT. 0) THEN
!**** Zonal Flow initialization *******************************************
! Every mode is zero
phi = 0._dp
moments_e = 0._dp
moments_i = 0._dp
! Except at ikr = mode number + 1, symmetry is already included since kr>=0
IF( (INIT_ZF+1 .GT. ikrs) .AND. (INIT_ZF+1 .LT. ikre) ) THEN
phi(INIT_ZF+1,ikz_0) = (2._dp*PI)**2*init_background/deltakr/deltakz/2._dp
ENDIF
!**** Init ZF in density
moments_e(1,1,INIT_ZF+1,ikz_0,:) = -(2._dp*PI)**2*init_background/deltakr/deltakz * imagu/2._dp
moments_i(1,1,INIT_ZF+1,ikz_0,:) = -(2._dp*PI)**2*init_background/deltakr/deltakz * imagu/2._dp
ELSE ! we compute phi from noisy moments and poisson
CALL poisson
ENDIF
END SUBROUTINE init_phi
!******************************************************************************!
!******************************************************************************!
!!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
!******************************************************************************!
SUBROUTINE build_dnjs_table
USE basic
USE array, Only : dnjs
USE grid, Only : jmaxe, jmaxi
USE coeff
IMPLICIT NONE
INTEGER :: in, ij, is, J
INTEGER :: n_, j_, s_
J = max(jmaxe,jmaxi)
DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
n_ = in - 1
DO ij = in,J+1
j_ = ij - 1
DO is = ij,J+1
s_ = is - 1
dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
! By symmetry
dnjs(in,is,ij) = dnjs(in,ij,is)
dnjs(ij,in,is) = dnjs(in,ij,is)
dnjs(ij,is,in) = dnjs(in,ij,is)
dnjs(is,ij,in) = dnjs(in,ij,is)
dnjs(is,in,ij) = dnjs(in,ij,is)
ENDDO
ENDDO
ENDDO
END SUBROUTINE build_dnjs_table
!******************************************************************************!
!******************************************************************************!
!!!!!!! Evaluate the kernels once for all
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
USE array, Only : kernel_e, kernel_i
USE grid
use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS
IMPLICIT NONE
REAL(dp) :: factj, j_dp, j_int
REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2
REAL(dp) :: be_2, bi_2, alphaD
REAL(dp) :: kr, kz, kperp2
!!!!! Electron kernels !!!!!
!Precompute species dependant factors
sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
factj = 1.0 ! Start of the recursive factorial
DO ij = 1, jmaxe+1
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp) ! REAL of degree
! Recursive factorial
IF (j_dp .GT. 0) THEN
factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
DO ikr = ikrs,ikre
kr = krarray(ikr)
DO ikz = ikzs,ikze
kz = kzarray(ikz)
be_2 = (kr**2 + kz**2) * sigmae2_taue_o2
kernel_e(ij, ikr, ikz) = be_2**j_int * exp(-be_2)/factj
ENDDO
ENDDO
ENDDO
! Kernels closure
DO ikr = ikrs,ikre
kr = krarray(ikr)
DO ikz = ikzs,ikze
kz = kzarray(ikz)
be_2 = (kr**2 + kz**2) * sigmae2_taue_o2
kernel_e(ijeg_e,ikr,ikz) = be_2/(real(ijeg_e,dp))*kernel_e(ije_e,ikr,ikz)
ENDDO
ENDDO
!!!!! Ion kernels !!!!!
sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
factj = 1.0 ! Start of the recursive factorial
DO ij = 1, jmaxi+1
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp) ! REAL of degree
! Recursive factorial
IF (j_dp .GT. 0) THEN
factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
DO ikr = ikrs,ikre
kr = krarray(ikr)
DO ikz = ikzs,ikze
kz = kzarray(ikz)
bi_2 = (kr**2 + kz**2) * sigmai2_taui_o2
kernel_i(ij, ikr, ikz) = bi_2**j_int * exp(-bi_2)/factj
ENDDO
ENDDO
ENDDO
! Kernels closure
DO ikr = ikrs,ikre
kr = krarray(ikr)
DO ikz = ikzs,ikze
kz = kzarray(ikz)
bi_2 = (kr**2 + kz**2) * sigmai2_taui_o2
kernel_i(ijeg_i,ikr,ikz) = bi_2/(real(ijeg_i,dp))*kernel_e(ije_i,ikr,ikz)
ENDDO
ENDDO
END SUBROUTINE evaluate_kernels
!******************************************************************************!
c4science · Help