-
Antoine Cyril David Hoffmann authoredAntoine Cyril David Hoffmann authored
inital.F90 10.33 KiB
!******************************************************************************!
!!!!!! 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 time_integration
USE array, ONLY : Sepj,Sipj
USE collision
USE closure
USE ghosts
USE restarts
USE numerics, ONLY: wipe_turbulence, wipe_zonalflow
IMPLICIT NONE
CALL set_updatetlevel(1)
!!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
! through loading a previous state
IF ( job2load .GE. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
CALL load_moments ! get N_0
CALL poisson ! compute phi_0=phi(N_0)
! through initialization
ELSE
! set phi with noise and set moments to 0
IF (INIT_NOISY_PHI) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
! set moments_00 (GC density) with noise and compute phi afterwards
ELSE
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
CALL init_gyrodens ! init only gyrocenter density
! CALL init_moments ! init all moments randomly (unadvised)
CALL poisson ! get phi_0 = phi(N_0)
ENDIF
ENDIF
! Option for wiping the ZF modes (ky==0)
IF ( WIPE_ZF .GT. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) '-Wiping ZF modes'
CALL wipe_zonalflow
ENDIF
! Option for wiping the turbulence and check growth of secondary inst.
IF ( WIPE_TURB .GT. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) '-Wiping turbulence'
CALL wipe_turbulence
ENDIF
! Option for initializing a gaussian blob on the zonal profile
IF ( INIT_BLOB ) THEN
IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
CALL initialize_blob
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(*,*) '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 all 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) :: kx, ky, sigma, gain, ky_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 ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
END DO
IF ( contains_kx0 ) THEN
DO iky=2,Nky/2 !symmetry at kx = 0 for all z
moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
END DO
ENDIF
END DO
END DO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
END DO
IF ( contains_kx0 ) THEN
DO iky=2,Nky/2 !symmetry at kx = 0 for all z
moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
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 ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE init_moments
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize the gyrocenter density randomly
!******************************************************************************!
SUBROUTINE init_gyrodens
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) :: kx, ky, sigma, gain, ky_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 ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE
moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
ENDIF
END DO
END DO
END DO
IF ( contains_kx0 ) THEN
DO iky=2,Nky/2 !symmetry at kx = 0 for all z
moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
END DO
ENDIF
END DO
END DO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE
moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
ENDIF
END DO
END DO
END DO
IF ( contains_kx0 ) THEN
DO iky=2,Nky/2 !symmetry at kx = 0 for all z
moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
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 ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE init_gyrodens
!******************************************************************************!
!******************************************************************************!
!!!!!!! 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) :: kx, ky, sigma, gain, ky_shift
INTEGER, DIMENSION(12) :: iseedarr
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr+my_id)
!**** noise initialization *******************************************
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
ENDDO
END DO
END DO
!symmetry at kx = 0 to keep real inverse transform
IF ( contains_kx0 ) THEN
DO iky=2,Nky/2
phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
END DO
phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real
ENDIF
!**** ensure no previous moments initialization
moments_e = 0._dp; moments_i = 0._dp
!**** Zonal Flow initialization *******************************************
! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
IF(INIT_ZF .GT. 0) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
DO iz = izs,ize
phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
ENDDO
ENDIF
ENDIF
END SUBROUTINE init_phi
!******************************************************************************!
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
!******************************************************************************!
SUBROUTINE initialize_blob
USE fields
USE grid
USE model, ONLY: sigmai2_taui_o2
IMPLICIT NONE
REAL(dp) ::kx, ky, sigma, gain
sigma = 0.5_dp
gain = 5e2_dp
DO ikx=ikxs,ikxe
kx = kxarray(ikx)
DO iky=ikys,ikye
ky = kyarray(iky)
DO iz=izs,ize
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :) &
+ gain*sigma/SQRT2 * exp(-(kx**2+ky**2)*sigma**2/4._dp) &
* AA_x(ikx)*AA_y(iky)!&
! * exp(sigmai2_taui_o2*(kx**2+ky**2))
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
END SUBROUTINE initialize_blob
!******************************************************************************!