From faea4734d906df3aff8d8e1427a9c5b531ed1f4a Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 26 May 2023 15:35:18 +0200
Subject: [PATCH] surprise
---
src/inital.F90 | 88 +++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 87 insertions(+), 1 deletion(-)
diff --git a/src/inital.F90 b/src/inital.F90
index 6a1eca0e..3dd50a2a 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -66,7 +66,13 @@ SUBROUTINE inital
CALL update_ghosts_moments
CALL solve_EM_fields
CALL update_ghosts_EM
- END SELECT
+ CASE('ricci')
+ CALL speak('Init Ricci')
+ CALL init_ricci ! init only gyrocenter density
+ 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
CALL speak('Apply closure')
@@ -448,3 +454,83 @@ SUBROUTINE init_ppj
ENDDO
END SUBROUTINE init_ppj
!******************************************************************************!
+
+!******************************************************************************!
+!!!!!!! Initialize Ricci density
+!******************************************************************************!
+SUBROUTINE init_ricci
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ local_nkx_offset, local_nky_offset,&
+ ngp, ngj, ngz, iky0, parray, jarray, contains_ky0, AA_x, AA_y
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp, imagu
+ USE initial_par,ONLY: iseed, init_noiselvl, init_background
+ USE model, ONLY: LINEARITY
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+
+ COMPLEX(xp), DIMENSION(186,52) :: ricci_mat_real, ricci_mat_imag, ricci_face
+ REAL(xp) :: tmp_real, tmp_imag
+ COMPLEX(xp) :: scaling
+ INTEGER :: ia,ip,ij,ikx,iky,iz
+ ! open data file
+ ricci_mat_real = 0; ricci_mat_imag = 0
+ open(unit = 1 , file = "/home/ahoffman/gyacomo/Gallery/fourier_ricci_real.txt")
+ read(1,*) ricci_mat_real
+ close(1)
+ open(unit = 2 , file = "/home/ahoffman/gyacomo/Gallery/fourier_ricci_imag.txt")
+ read(2,*) ricci_mat_imag
+ close(2)
+ scaling = 0.000002
+ moments = 0._xp
+ !**** Broad noise initialization *******************************************
+ DO ia=1,local_na
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO ip=1+ngp/2,local_np+ngp/2
+ DO ij=1+ngj/2,local_nj+ngj/2
+ DO iz=1+ngz/2,local_nz+ngz/2
+ IF((ikx+local_nkx_offset .LE. 186) .AND. (iky+local_nky_offset .LE. 52)) THEN
+ IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
+ moments(ia,ip,ij,iky,ikx,iz,:) = scaling*(ricci_mat_real(ikx+local_nkx_offset,iky+local_nky_offset)&
+ + imagu*ricci_mat_imag(ikx+local_nkx_offset,iky+local_nky_offset))
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ENDIF
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ENDIF
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ print*, sum(moments)
+ IF ( contains_ky0 ) THEN
+ DO ip=1+ngp/2,local_np+ngp/2
+ DO ij=1+ngj/2,local_nj+ngj/2
+ DO iz=1+ngz/2,local_nz+ngz/2
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
+ moments(ia, ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:,:)
+ END DO
+ END DO
+ END DO
+ END DO
+ ENDIF
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+END SUBROUTINE init_ricci
+!******************************************************************************!
\ No newline at end of file
--
GitLab