Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • dev
  • ExB_shearing
  • one_moment_array
  • gyacomo22
  • ahoffman-master-patch-08280
  • Miller
  • shear_version
  • gather_output
  • divided_output
  • half_ky_plane
  • half_kx_plane
  • 3D_parallel
  • HeLaZ2
  • 2D
  • 3D
  • PPB110
  • v3.0.3
  • V3.0.2
  • V2.1.0
  • v3.0.1
  • v2.0
  • v3.0
23 results
Blame Permalink
inital.F90 16.20 KiB 
!******************************************************************************!
!!!!!! initialize the moments and load/build coeff tables
!******************************************************************************!
SUBROUTINE inital

  USE basic,            ONLY: my_id
  USE initial_par,      ONLY: INIT_OPT
  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_EM
  USE restarts,         ONLY: load_moments, job2load
  USE processing,       ONLY: compute_fluid_moments
  USE model,            ONLY: KIN_E, LINEARITY
  USE nonlinear,        ONLY: compute_Sapj, nonlinear_init
  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 update_ghosts_moments
    CALL solve_EM_fields ! compute phi_0=phi(N_0)
    CALL update_ghosts_EM
  ! through initialization
  ELSE
    SELECT CASE (INIT_OPT)
    ! set phi with noise and set moments to 0
    CASE ('phi')
      IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
      CALL init_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 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 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 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 solve_EM_fields
      CALL update_ghosts_EM
    END SELECT
  ENDIF
  ! closure of j>J, p>P and j<0, p<0 moments
  IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
  CALL apply_closure_model
  ! ghosts for p parallelization
  IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
  CALL update_ghosts_moments
  CALL update_ghosts_EM
  !! End of phi and moments initialization

  ! Load the COSOlver collision operator coefficients
  IF(cosolver_coll) &
  CALL load_COSOlver_mat

  !! Preparing auxiliary arrays at initial state
  ! particle density, fluid velocity and temperature (used in diagnose)
  IF (my_id .EQ. 0) WRITE(*,*) 'Computing fluid moments'
  CALL compute_fluid_moments

  ! init auxval for nonlinear
  CALL nonlinear_init
  ! compute nonlinear for t=0 diagnostic
  CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)

END SUBROUTINE inital
!******************************************************************************!

!******************************************************************************!
!!!!!!! Initialize all the moments randomly
!******************************************************************************!
SUBROUTINE init_moments
  USE basic
  USE grid
  USE initial_par,ONLY: iseed, init_noiselvl, init_background
  USE fields,     ONLY: moments_e, moments_i
  USE prec_const, ONLY: dp
  USE utility,    ONLY: checkfield
  USE model,      ONLY : LINEARITY, KIN_E
  IMPLICIT NONE

  REAL(dp) :: noise
  INTEGER, DIMENSION(12) :: iseedarr

  ! Seed random number generator
  iseedarr(:)=iseed
  CALL RANDOM_SEED(PUT=iseedarr+my_id)

    !**** Broad noise initialization *******************************************
    ! Electron init
    IF(KIN_E) THEN
    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,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
            END DO
          END DO
        END DO
        IF ( contains_ky0 ) THEN
          DO iky=2,Nky/2 !symmetry at ky = 0 for all z
            moments_e(ip,ij,iky_0,ikx,:,:) = moments_e( ip,ij,iky_0,Nkx+2-ikx,:, :)
          END DO
        ENDIF
      END DO
    END DO
    ENDIF
    ! Ion init
    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,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
            END DO
          END DO
        END DO

        IF ( contains_ky0 ) THEN
          DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
            moments_i( ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:,:)
          END DO
        ENDIF

      END DO
    END DO

    ! Putting to zero modes that are not in the 2/3 Orszag rule
    IF (LINEARITY .NE. 'linear') THEN
      DO ikx=ikxs,ikxe
      DO iky=ikys,ikye
      DO iz=izs,ize
        IF(KIN_E) THEN
        DO ip=ips_e,ipe_e
        DO ij=ijs_e,ije_e
          moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
        ENDDO
        ENDDO
        ENDIF
        DO ip=ips_i,ipe_i
        DO ij=ijs_i,ije_i
          moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,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: KIN_E, LINEARITY
  IMPLICIT NONE

  REAL(dp) :: noise
  INTEGER, DIMENSION(12) :: iseedarr

  ! Seed random number generator
  iseedarr(:)=iseed
  CALL RANDOM_SEED(PUT=iseedarr+my_id)

    !**** Broad noise initialization *******************************************
    IF(KIN_E) THEN
    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,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
              ELSE
                moments_e(ip,ij,iky,ikx,iz,:) = 0._dp
              ENDIF
            END DO
          END DO
        END DO
        IF ( contains_ky0 ) THEN
          DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
            moments_e(ip,ij,iky_0,ikx,:,:) = moments_e(ip,ij,iky_0,Nkx+2-ikx,:,:)
          END DO
        ENDIF
      END DO
    END DO
    ENDIF

    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,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
              ELSE
                moments_i(ip,ij,iky,ikx,iz,:) = 0._dp
              ENDIF
            END DO
          END DO
        END DO
        IF ( contains_ky0 ) THEN
          DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
            moments_i( ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:,:)
          END DO
        ENDIF
      END DO
    END DO

    ! Putting to zero modes that are not in the 2/3 Orszag rule
    IF (LINEARITY .NE. 'linear') THEN
      DO ikx=ikxs,ikxe
      DO iky=ikys,ikye
      DO iz=izs,ize
        IF(KIN_E) THEN
        DO ip=ips_e,ipe_e
        DO ij=ijs_e,ije_e
          moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
        ENDDO
        ENDDO
        ENDIF
        DO ip=ips_i,ipe_i
        DO ij=ijs_i,ije_i
          moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,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
  USE model, ONLY: KIN_E, LINEARITY
  IMPLICIT NONE

  REAL(dp) :: noise
  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(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
        ENDDO
      END DO
    END DO

    !symmetry at ky = 0 to keep real inverse transform
    IF ( contains_ky0 ) THEN
      DO ikx=2,Nkx/2
        phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize)
      END DO
      phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize)) !origin must be real
    ENDIF

    !**** ensure no previous moments initialization
    IF(KIN_E) 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(iky_0,INIT_ZF+1,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
          moments_i(1,1,iky_0,INIT_ZF+1,iz,:) = kxarray(INIT_ZF+1)**2*phi(iky_0,INIT_ZF+1,iz)* COS((iz-1)/Nz*2._dp*PI)
          IF(KIN_E) moments_e(1,1,iky_0,INIT_ZF+1,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: KIN_E, LINEARITY
  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,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,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
    IF(KIN_E) THEN
    DO ip=ips_e,ipe_e
    DO ij=ijs_e,ije_e
      IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
        moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,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
    ENDIF
  ENDDO
  ENDDO
  ENDDO
END SUBROUTINE initialize_blob
!******************************************************************************!


!******************************************************************************!
!!!!!!! Initialize the gyrocenter in a ppj gene manner (power law)
!******************************************************************************!
SUBROUTINE init_ppj
  USE basic
  USE grid
  USE fields, ONLY: moments_e, moments_i
  USE array,  ONLY: kernel_e, kernel_i
  USE prec_const
  USE utility, ONLY: checkfield
  USE initial_par
  USE model, ONLY: KIN_E, LINEARITY
  USE geometry, ONLY: Jacobian, iInt_Jacobian

  IMPLICIT NONE

  REAL(dp) :: kx, ky, sigma_z, amp, z

  sigma_z = pi/4._dp
  amp = 1.0_dp

    !**** Broad noise initialization *******************************************
    ! Electrons
    IF (KIN_E) THEN
    DO ip=ips_e,ipe_e
      DO ij=ijs_e,ije_e
        IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
          DO ikx=ikxs,ikxe
            kx = kxarray(ikx)
            DO iky=ikys,ikye
              ky = kyarray(iky)
              DO iz=izs,ize
                z = zarray(iz,0)
                IF (kx .EQ. 0) THEN
                  IF(ky .EQ. 0) THEN
                    moments_e(ip,ij,iky,ikx,iz,:) = 0._dp
                  ELSE
                    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,:) = (deltakx/(ABS(kx)+deltakx))*(ky_min/(ABS(ky)+ky_min))
                  ELSE
                    moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(deltakx/(ABS(kx)+deltakx))
                  ENDIF
                ENDIF
                ! z-dep and noise
                moments_e(ip,ij,iky,ikx,iz,:) = moments_e(ip,ij,iky,ikx,iz,:) * &
                (Jacobian(iz,0)*iInt_Jacobian)**2

                ! divide by kernel_0 to adjust to particle density (n = kernel_0 N00)
                moments_e(ip,ij,iky,ikx,iz,:) = moments_e(ip,ij,iky,ikx,iz,:)/kernel_e(ij,iky,ikx,iz,0)
              END DO
            END DO
          END DO

          IF ( contains_ky0 ) THEN
            DO ikx=2,Nkx/2 !symmetry at kx = 0 for all z
              moments_e(ip,ij,iky_0,ikx,:,:) = moments_e( ip,ij,iky_0,Nkx+2-ikx,:, :)
            END DO
          ENDIF
        ELSE
          moments_e(ip,ij,:,:,:,:) = 0._dp
        ENDIF
      END DO
    END DO
    ENDIF

    ! Ions
    DO ip=ips_i,ipe_i
      DO ij=ijs_i,ije_i
        IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
          DO ikx=ikxs,ikxe
            kx = kxarray(ikx)
            DO iky=ikys,ikye
              ky = kyarray(iky)
              DO iz=izs,ize
                z = zarray(iz,0)
                IF (kx .EQ. 0) THEN
                  IF(ky .EQ. 0) THEN
                    moments_i(ip,ij,iky,ikx,iz,:) = 0._dp
                  ELSE
                    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,:) = (deltakx/(ABS(kx)+deltakx))*(ky_min/(ABS(ky)+ky_min))
                  ELSE
                    moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(deltakx/(ABS(kx)+deltakx))
                  ENDIF
                ENDIF
                ! z-dep and noise
                moments_i(ip,ij,iky,ikx,iz,:) = moments_i(ip,ij,iky,ikx,iz,:) * &
                (Jacobian(iz,0)*iInt_Jacobian)**2
                ! divide by kernel_0 to adjust to particle density (n = kernel_0 N00)
                moments_i(ip,ij,iky,ikx,iz,:) = moments_i(ip,ij,iky,ikx,iz,:)/kernel_i(ij,iky,ikx,iz,0)
              END DO
            END DO
          END DO

          IF ( contains_ky0 ) THEN
            DO ikx=2,Nkx/2 !symmetry at kx = 0 for all z
              moments_i(ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:, :)
            END DO
          ENDIF
        ELSE
          moments_i(ip,ij,:,:,:,:) = 0._dp
        ENDIF
      END DO
    END DO

    ! Putting to zero modes that are not in the 2/3 Orszag rule
    IF (LINEARITY .NE. 'linear') 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,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,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,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
        ENDDO
        ENDDO
      ENDDO
      ENDDO
      ENDDO
    ENDIF

END SUBROUTINE init_ppj
!******************************************************************************!