From cbaf51c1d1e98ece931a7d62d157cf110e183203 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 31 Jan 2023 13:58:09 +0100
Subject: [PATCH] Improved blob init and added a phi ppj init
---
src/inital.F90 | 80 +++++++++++++++++++++++++++++++++++++++++++-------
1 file changed, 69 insertions(+), 11 deletions(-)
diff --git a/src/inital.F90 b/src/inital.F90
index c7c0247b..f7206588 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -33,6 +33,10 @@ SUBROUTINE inital
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
CALL update_ghosts_EM
+ CASE ('phi_ppj')
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
+ CALL init_phi_ppj
+ 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'
@@ -325,6 +329,54 @@ SUBROUTINE init_phi
END SUBROUTINE init_phi
!******************************************************************************!
+!******************************************************************************!
+!!!!!!! Initialize a ppj ES potential and cancel the moments
+!******************************************************************************!
+SUBROUTINE init_phi_ppj
+ USE basic
+ USE grid
+ USE fields
+ USE prec_const
+ USE initial_par
+ USE model, ONLY: KIN_E, LINEARITY
+ USE geometry, ONLY: Jacobian, iInt_Jacobian
+ IMPLICIT NONE
+ REAL(dp) :: kx, ky, z, amp
+ amp = 1.0_dp
+ !**** ppj initialization *******************************************
+ DO ikx=ikxs,ikxe
+ kx = kxarray(ikx)
+ DO iky=ikys,ikye
+ ky = kyarray(iky)
+ DO iz=izs,ize
+ z = zarray(iz,0)
+ IF (ky .NE. 0) THEN
+ phi(iky,ikx,iz) = 0._dp
+ ELSE
+ phi(iky,ikx,iz) = 0.5_dp*amp*(deltakx/(ABS(kx)+deltakx))
+ ENDIF
+ ! z-dep and noise
+ phi(iky,ikx,iz) = phi(iky,ikx,iz) * &
+ (Jacobian(iz,0)*iInt_Jacobian)**2
+ END DO
+ 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
+
+END SUBROUTINE init_phi_ppj
+!******************************************************************************!
+
!******************************************************************************!
!******************************************************************************!
!!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
@@ -332,23 +384,28 @@ END SUBROUTINE init_phi
SUBROUTINE initialize_blob
USE fields
USE grid
+ USE geometry, ONLY: shear, Jacobian, iInt_Jacobian
USE model, ONLY: KIN_E, LINEARITY
IMPLICIT NONE
- REAL(dp) ::kx, ky, sigma, gain
- sigma = 0.5_dp
- gain = 5e2_dp
+ REAL(dp) ::kx, ky, z, sigma_x, sigma_y, gain
+ sigma_y = 1.0
+ sigma_x = sigma_y
+ gain = 10.0
+
- DO ikx=ikxs,ikxe
- kx = kxarray(ikx)
DO iky=ikys,ikye
ky = kyarray(iky)
DO iz=izs,ize
+ z = zarray(iz,0)
+ DO ikx=ikxs,ikxe
+ kx = kxarray(ikx) + z*shear*ky
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)!&
+ + gain * exp(-((kx/sigma_x)**2+(ky/sigma_y)**2)) &
+ * AA_x(ikx)*AA_y(iky)* &
+ (Jacobian(iz,0)*iInt_Jacobian)**2!&
! * exp(sigmai2_taui_o2*(kx**2+ky**2))
ENDIF
ENDDO
@@ -358,8 +415,9 @@ SUBROUTINE initialize_blob
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)!&
+ + gain * exp(-((kx/sigma_x)**2+(ky/sigma_y)**2)) &
+ * AA_x(ikx)*AA_y(iky)* &
+ (Jacobian(iz,0)*iInt_Jacobian)**2!&
! * exp(sigmai2_taui_o2*(kx**2+ky**2))
ENDIF
ENDDO
@@ -423,7 +481,7 @@ SUBROUTINE init_ppj
(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)
+ ! 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
@@ -467,7 +525,7 @@ SUBROUTINE init_ppj
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)
+ ! 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
--
GitLab