diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90 index b182224..81d23ab 100644 --- a/src/ghosts_mod.F90 +++ b/src/ghosts_mod.F90 @@ -1,304 +1,312 @@ module ghosts USE basic USE fields, ONLY : moments_e, moments_i, phi USE grid USE time_integration USE model, ONLY : KIN_E USE geometry, ONLY : SHEARED, ikx_zBC_map IMPLICIT NONE INTEGER :: status(MPI_STATUS_SIZE), source, dest, count, ipg -PUBLIC :: update_ghosts +PUBLIC :: update_ghosts_moments, update_ghosts_phi CONTAINS -SUBROUTINE update_ghosts +SUBROUTINE update_ghosts_moments CALL cpu_time(t0_ghost) IF (num_procs_p .GT. 1) THEN ! Do it only if we share the p IF(KIN_E)& CALL update_ghosts_p_e CALL update_ghosts_p_i ENDIF IF(Nz .GT. 1) THEN IF(KIN_E) & CALL update_ghosts_z_e CALL update_ghosts_z_i - CALL update_ghosts_z_phi ENDIF tc_ghost = tc_ghost + (t1_ghost - t0_ghost) -END SUBROUTINE update_ghosts +END SUBROUTINE update_ghosts_moments + +SUBROUTINE update_ghosts_phi + CALL cpu_time(t0_ghost) + + IF(Nz .GT. 1) THEN + CALL update_ghosts_z_phi + ENDIF + tc_ghost = tc_ghost + (t1_ghost - t0_ghost) +END SUBROUTINE update_ghosts_phi !Communicate p+1, p+2 moments to left neighboor and p-1, p-2 moments to right one ! [a b|C D|e f] : proc n has moments a to f where a,b,e,f are ghosts ! !proc 0: [0 1 2 3 4|5 6] ! V V ^ ^ !proc 1: [3 4|5 6 7 8|9 10] ! V V ^ ^ !proc 2: [7 8|9 10 11 12|13 14] ! V V ^ ^ !proc 3: [11 12|13 14 15 16|17 18] ! ^ ^ !Closure by zero truncation : 0 0 SUBROUTINE update_ghosts_p_e IMPLICIT NONE count = (ijge_e-ijgs_e+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1) !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!! ! Send the last local moment to fill the -1 neighbour ghost CALL mpi_sendrecv(moments_e(ipe_e ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 10, & ! Send to right moments_e(ips_e-1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 10, & ! Recieve from left comm0, status, ierr) IF (deltape .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil CALL mpi_sendrecv(moments_e(ipe_e-1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 11, & ! Send to right moments_e(ips_e-2,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 11, & ! Recieve from left comm0, status, ierr) !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(moments_e(ips_e ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 12, & ! Send to left moments_e(ipe_e+1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 12, & ! Recieve from right comm0, status, ierr) IF (deltape .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil CALL mpi_sendrecv(moments_e(ips_e+1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 13, & ! Send to left moments_e(ipe_e+2,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 13, & ! Recieve from right comm0, status, ierr) END SUBROUTINE update_ghosts_p_e !Communicate p+1, p+2 moments to left neighboor and p-1, p-2 moments to right one SUBROUTINE update_ghosts_p_i IMPLICIT NONE count = (ijge_i-ijgs_i+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1) ! Number of elements sent !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(moments_i(ipe_i ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 14, & moments_i(ips_i-1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 14, & comm0, status, ierr) IF (deltapi .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil CALL mpi_sendrecv(moments_i(ipe_i-1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 15, & moments_i(ips_i-2,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 15, & comm0, status, ierr) !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(moments_i(ips_i ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 16, & moments_i(ipe_i+1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 16, & comm0, status, ierr) IF (deltapi .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil CALL mpi_sendrecv(moments_i(ips_i+1,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 17, & moments_i(ipe_i+2,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 17, & comm0, status, ierr) END SUBROUTINE update_ghosts_p_i !Communicate z+1, z+2 moments to left neighboor and z-1, z-2 moments to right one ! [a b|C D|e f] : proc n has moments a to f where a,b,e,f are ghosts ! !proc 0: [0 1 2 3 4|5 6] ! V V ^ ^ !proc 1: [3 4|5 6 7 8|9 10] ! V V ^ ^ !proc 2: [7 8|9 10 11 12|13 14] ! V V ^ ^ !proc 3: [11 12|13 14 15 16|17 18] ! ^ ^ !Periodic: 0 1 SUBROUTINE update_ghosts_z_e IMPLICIT NONE INTEGER :: ikxBC CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) IF (num_procs_z .GT. 1) THEN count = (ipge_e-ipgs_e+1)*(ijge_e-ijgs_e+1)*(ikye-ikys+1)*(ikxe-ikxs+1) !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! ! Send the last local moment to fill the -1 neighbour ghost CALL mpi_sendrecv(moments_e(:,:,:,:,ize ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 20, & ! Send to Up the last moments_e(:,:,:,:,izs-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 20, & ! Recieve from Down the first-1 comm0, status, ierr) CALL mpi_sendrecv(moments_e(:,:,:,:,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 21, & ! Send to Up the last-1 moments_e(:,:,:,:,izs-2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 21, & ! Recieve from Down the first-2 comm0, status, ierr) !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(moments_e(:,:,:,:,izs ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 22, & ! Send to Down the first moments_e(:,:,:,:,ize+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 22, & ! Recieve from Up the last+1 comm0, status, ierr) CALL mpi_sendrecv(moments_e(:,:,:,:,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 23, & ! Send to Down the first+1 moments_e(:,:,:,:,ize+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 23, & ! Recieve from Up the last+2 comm0, status, ierr) ELSE ! still need to perform periodic boundary conditions IF(SHEARED) THEN DO iky = ikys,ikye DO ikx = ikxs,ikxe ikxBC = ikx_zBC_map(ikx,iky); IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a) ! first-1 gets last moments_e(iky,ikx,:,:,izs-1,updatetlevel) = moments_e(iky,ikxBC,:,:,ize ,updatetlevel) ! first-2 gets last-1 moments_e(iky,ikx,:,:,izs-2,updatetlevel) = moments_e(iky,ikxBC,:,:,ize-1,updatetlevel) ! last+1 gets first moments_e(iky,ikx,:,:,ize+1,updatetlevel) = moments_e(iky,ikxBC,:,:,izs ,updatetlevel) ! last+2 gets first+1 moments_e(iky,ikx,:,:,ize+2,updatetlevel) = moments_e(iky,ikxBC,:,:,izs+1,updatetlevel) ELSE moments_e(iky,ikx,:,:,izs-1,updatetlevel) = 0._dp moments_e(iky,ikx,:,:,izs-2,updatetlevel) = 0._dp moments_e(iky,ikx,:,:,ize+1,updatetlevel) = 0._dp moments_e(iky,ikx,:,:,ize+2,updatetlevel) = 0._dp ENDIF ENDDO ENDDO ELSE ! No shear so simple periodic BC ! first-1 gets last moments_e(:,:,:,:,izs-1,updatetlevel) = moments_e(:,:,:,:,ize ,updatetlevel) ! first-2 gets last-1 moments_e(:,:,:,:,izs-2,updatetlevel) = moments_e(:,:,:,:,ize-1,updatetlevel) ! last+1 gets first moments_e(:,:,:,:,ize+1,updatetlevel) = moments_e(:,:,:,:,izs ,updatetlevel) ! last+2 gets first+1 moments_e(:,:,:,:,ize+2,updatetlevel) = moments_e(:,:,:,:,izs+1,updatetlevel) ENDIF ENDIF END SUBROUTINE update_ghosts_z_e SUBROUTINE update_ghosts_z_i IMPLICIT NONE INTEGER :: ikxBC CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) IF (num_procs_z .GT. 1) THEN count = (ipge_i-ipgs_i+1)*(ijge_i-ijgs_i+1)*(ikye-ikys+1)*(ikxe-ikxs+1) !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! ! Send the last local moment to fill the -1 neighbour ghost CALL mpi_sendrecv(moments_i(:,:,:,:,ize ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 30, & ! Send to Up the last moments_i(:,:,:,:,izs-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 30, & ! Recieve from Down the first-1 comm0, status, ierr) CALL mpi_sendrecv(moments_i(:,:,:,:,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 31, & ! Send to Up the last-1 moments_i(:,:,:,:,izs-2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 31, & ! Recieve from Down the first-2 comm0, status, ierr) !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(moments_i(:,:,:,:,izs ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 32, & ! Send to Down the first moments_i(:,:,:,:,ize+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 32, & ! Recieve from Down the last+1 comm0, status, ierr) CALL mpi_sendrecv(moments_i(:,:,:,:,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 33, & ! Send to Down the first+1 moments_i(:,:,:,:,ize+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 33, & ! Recieve from Down the last+2 comm0, status, ierr) ELSE ! still need to perform periodic boundary conditions IF(SHEARED) THEN DO iky = ikys,ikye DO ikx = ikxs,ikxe ikxBC = ikx_zBC_map(iky,ikx); IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a) ! first-1 gets last moments_i(:,:,iky,ikx,izs-1,updatetlevel) = moments_i(:,:,iky,ikxBC,ize ,updatetlevel) ! first-2 gets last-1 moments_i(:,:,iky,ikx,izs-2,updatetlevel) = moments_i(:,:,iky,ikxBC,ize-1,updatetlevel) ! last+1 gets first moments_i(:,:,iky,ikx,ize+1,updatetlevel) = moments_i(:,:,iky,ikxBC,izs ,updatetlevel) ! last+2 gets first+1 moments_i(:,:,iky,ikx,ize+2,updatetlevel) = moments_i(:,:,iky,ikxBC,izs+1,updatetlevel) ELSE moments_i(:,:,iky,ikx,izs-1,updatetlevel) = 0._dp moments_i(:,:,iky,ikx,izs-2,updatetlevel) = 0._dp moments_i(:,:,iky,ikx,ize+1,updatetlevel) = 0._dp moments_i(:,:,iky,ikx,ize+2,updatetlevel) = 0._dp ENDIF ENDDO ENDDO ELSE ! No shear so simple periodic BC ! first-1 gets last moments_i(:,:,:,:,izs-1,updatetlevel) = moments_i(:,:,:,:,ize ,updatetlevel) ! first-2 gets last-1 moments_i(:,:,:,:,izs-2,updatetlevel) = moments_i(:,:,:,:,ize-1,updatetlevel) ! last+1 gets first moments_i(:,:,:,:,ize+1,updatetlevel) = moments_i(:,:,:,:,izs ,updatetlevel) ! last+2 gets first+1 moments_i(:,:,:,:,ize+2,updatetlevel) = moments_i(:,:,:,:,izs+1,updatetlevel) ENDIF ENDIF END SUBROUTINE update_ghosts_z_i SUBROUTINE update_ghosts_z_phi IMPLICIT NONE INTEGER :: ikxBC CALL cpu_time(t1_ghost) IF(Nz .GT. 1) THEN CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) IF (num_procs_z .GT. 1) THEN count = (ikye-ikys+1) * (ikxe-ikxs+1) !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! ! Send the last local moment to fill the -1 neighbour ghost CALL mpi_sendrecv(phi(:,:,ize ), count, MPI_DOUBLE_COMPLEX, nbr_U, 40, & ! Send to Up the last phi(:,:,izs-1), count, MPI_DOUBLE_COMPLEX, nbr_D, 40, & ! Receive from Down the first-1 comm0, status, ierr) CALL mpi_sendrecv(phi(:,:,ize-1), count, MPI_DOUBLE_COMPLEX, nbr_U, 41, & ! Send to Up the last-1 phi(:,:,izs-2), count, MPI_DOUBLE_COMPLEX, nbr_D, 41, & ! Receive from Down the first-2 comm0, status, ierr) !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(phi(:,:,izs ), count, MPI_DOUBLE_COMPLEX, nbr_D, 42, & ! Send to Down the first phi(:,:,ize+1), count, MPI_DOUBLE_COMPLEX, nbr_U, 42, & ! Recieve from Up the last+1 comm0, status, ierr) CALL mpi_sendrecv(phi(:,:,izs+1), count, MPI_DOUBLE_COMPLEX, nbr_D, 43, & ! Send to Down the first+1 phi(:,:,ize+2), count, MPI_DOUBLE_COMPLEX, nbr_U, 43, & ! Recieve from Up the last+2 comm0, status, ierr) ELSE ! still need to perform periodic boundary conditions phi(:,:,izs-1) = phi(:,:,ize) phi(:,:,izs-2) = phi(:,:,ize-1) phi(:,:,ize+1) = phi(:,:,izs) phi(:,:,ize+2) = phi(:,:,izs+1) IF(SHEARED) THEN DO iky = ikys,ikye DO ikx = ikxs,ikxe ikxBC = ikx_zBC_map(iky,ikx); IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a) ! first-1 gets last phi(iky,ikx,izs-1) = phi(iky,ikxBC,ize ) ! first-2 gets last-1 phi(iky,ikx,izs-2) = phi(iky,ikxBC,ize-1) ! last+1 gets first phi(iky,ikx,ize+1) = phi(iky,ikxBC,izs ) ! last+2 gets first+1 phi(iky,ikx,ize+2) = phi(iky,ikxBC,izs+1) ELSE phi(iky,ikx,izs-1) = 0._dp phi(iky,ikx,izs-2) = 0._dp phi(iky,ikx,ize+1) = 0._dp phi(iky,ikx,ize+2) = 0._dp ENDIF ENDDO ENDDO ELSE ! No shear so simple periodic BC ! first-1 gets last phi(:,:,izs-1) = phi(:,:,ize ) ! first-2 gets last-1 phi(:,:,izs-2) = phi(:,:,ize-1) ! last+1 gets first phi(:,:,ize+1) = phi(:,:,izs ) ! last+2 gets first+1 phi(:,:,ize+2) = phi(:,:,izs+1) ENDIF ENDIF ENDIF CALL cpu_time(t1_ghost) tc_ghost = tc_ghost + (t1_ghost - t0_ghost) END SUBROUTINE update_ghosts_z_phi END MODULE ghosts diff --git a/src/inital.F90 b/src/inital.F90 index 343c4f6..c692fb1 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,511 +1,517 @@ !******************************************************************************! !!!!!! 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 + USE ghosts, ONLY: update_ghosts_moments, update_ghosts_phi USE restarts, ONLY: load_moments, job2load USE numerics, ONLY: play_with_modes, save_EM_ZF_modes 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 poisson ! compute phi_0=phi(N_0) - CALL update_ghosts + CALL update_ghosts_phi ! 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 + CALL update_ghosts_phi ! 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 + CALL update_ghosts_moments CALL poisson + CALL update_ghosts_phi ! 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 + CALL update_ghosts_moments CALL poisson + CALL update_ghosts_phi ! init a gaussian blob in gyrodens CASE('blob') IF (my_id .EQ. 0) WRITE(*,*) '--init a blob' CALL initialize_blob - CALL update_ghosts + CALL update_ghosts_moments CALL poisson + CALL update_ghosts_phi ! 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 + CALL update_ghosts_moments CALL poisson + CALL update_ghosts_phi 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 + CALL update_ghosts_moments + CALL update_ghosts_phi !! End of phi and moments initialization ! Save (kx,0) and (0,ky) modes for num exp CALL save_EM_ZF_modes ! Freeze/Wipe some selected modes (entropy,zonal,turbulent) CALL play_with_modes ! 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 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 ******************************************* 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 iky=2,Nky/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 REAL(dp) :: kx, ky, kp, 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(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: sigmai2_taui_o2, 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 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) :: noise REAL(dp) :: kx, ky, sigma_z, amp, ky_shift, z INTEGER, DIMENSION(12) :: iseedarr sigma_z = pi/4.0 amp = 1e4 !**** 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._dp!0.5_dp * ky_min/(ABS(ky)+ky_min) ENDIF ELSE IF(ky .GT. 0) THEN moments_e(ip,ij,iky,ikx,iz,:) = 0._dp!(kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min)) ELSE moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(kx_min/(ABS(kx)+kx_min)) ENDIF ENDIF CALL RANDOM_NUMBER(noise) ! z-dep and noise moments_e(ip,ij,iky,ikx,iz,:) = moments_e(ip,ij,iky,ikx,iz,:) * & (Jacobian(iz,0)*iInt_Jacobian)**2 !+ (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 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._dp!0.5_dp * ky_min/(ABS(ky)+ky_min) ENDIF ELSE IF(ky .GT. 0) THEN moments_i(ip,ij,iky,ikx,iz,:) = 0._dp!(kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min)) ELSE moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp*amp*(kx_min/(ABS(kx)+kx_min)) ENDIF ENDIF CALL RANDOM_NUMBER(noise) ! z-dep and noise moments_i(ip,ij,iky,ikx,iz,:) = moments_i(ip,ij,iky,ikx,iz,:) * & (Jacobian(iz,0)*iInt_Jacobian)**2 !+ (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 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 ! Adjust the scaling to trigger faster NL saturation IF(KIN_E) & moments_e = 1e3*moments_e moments_i = 1e3*moments_i END SUBROUTINE init_ppj !******************************************************************************! diff --git a/src/stepon.F90 b/src/stepon.F90 index 7c4ee5f..0974971 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,165 +1,171 @@ SUBROUTINE stepon ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme) USE advance_field_routine, ONLY: advance_time_level, advance_field, advance_moments USE basic USE closure - USE ghosts, ONLY: update_ghosts + USE ghosts, ONLY: update_ghosts_moments, update_ghosts_phi USE grid USE model, ONLY : LINEARITY, KIN_E use prec_const USE time_integration USE numerics, ONLY: play_with_modes USE utility, ONLY: checkfield IMPLICIT NONE INTEGER :: num_step LOGICAL :: mlend DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4 !----- BEFORE: All fields+ghosts are updated for step = n ! Compute right hand side from current fields ! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n) CALL assemble_RHS ! ---- step n -> n+1 transition ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme) CALL advance_time_level ! ---- ! Update moments with the hierarchy RHS (step by step) ! N_n+1 = N_n + N_rhs(n) CALL advance_moments ! Closure enforcement of moments CALL apply_closure_model ! Exchanges the ghosts values of N_n+1 - CALL update_ghosts + CALL update_ghosts_moments ! Update electrostatic potential phi_n = phi(N_n+1) CALL poisson - CALL update_ghosts_z_phi + CALL update_ghosts_phi ! Numerical experiments ! Store or cancel/maintain zonal modes artificially ! CALL play_with_modes !- Check before next step CALL checkfield_all() IF( nlend ) EXIT ! exit do loop CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) !----- AFTER: All fields are updated for step = n+1 END DO CONTAINS !!!! Basic structure to simplify stepon SUBROUTINE assemble_RHS USE moments_eq_rhs, ONLY: compute_moments_eq_rhs USE collision, ONLY: compute_TColl USE nonlinear, ONLY: compute_Sapj USE processing, ONLY: compute_nadiab_moments_z_gradients_and_interp IMPLICIT NONE ! compute auxiliary non adiabatic moments array, gradients and interp CALL compute_nadiab_moments_z_gradients_and_interp ! compute nonlinear term ("if linear" is included inside) CALL compute_Sapj ! compute collision term ("if coll, if nu >0" is included inside) CALL compute_TColl ! compute the moments equation rhs CALL compute_moments_eq_rhs END SUBROUTINE assemble_RHS SUBROUTINE checkfield_all ! Check all the fields for inf or nan + USE fields, ONLY: phi, moments_e, moments_i + IMPLICIT NONE ! Execution time start CALL cpu_time(t0_checkfield) IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0 mlend=.FALSE. IF(.NOT.nlend) THEN mlend=mlend .or. checkfield(phi,' phi') IF(KIN_E) THEN DO ij=ijgs_e,ijge_e DO ip=ipgs_e,ipge_e mlend=mlend .or. checkfield(moments_e(ip,ij,:,:,:,updatetlevel),' moments_e') ENDDO ENDDO ENDIF DO ij=ijgs_i,ijge_i DO ip=ipgs_i,ipge_i mlend=mlend .or. checkfield(moments_i(ip,ij,:,:,:,updatetlevel),' moments_i') ENDDO ENDDO CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) ENDIF ! Execution time end CALL cpu_time(t1_checkfield) tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield) END SUBROUTINE checkfield_all SUBROUTINE anti_aliasing + USE fields, ONLY: moments_e, moments_i + IMPLICIT NONE IF(KIN_E)THEN DO iz=izgs,izge DO ikx=ikxs,ikxe DO iky=ikys,ikye DO ij=ijgs_e,ijge_e DO ip=ipgs_e,ipge_e moments_e( ip,ij,iky,ikx,iz,:) = AA_x(ikx)* AA_y(iky) * moments_e( ip,ij,iky,ikx,iz,:) END DO END DO END DO END DO END DO ENDIF DO iz=izgs,izge DO ikx=ikxs,ikxe DO iky=ikys,ikye DO ij=ijs_i,ije_i DO ip=ipgs_i,ipge_i moments_i( ip,ij,iky,ikx,iz,:) = AA_x(ikx)* AA_y(iky) * moments_i( ip,ij,iky,ikx,iz,:) END DO END DO END DO END DO END DO END SUBROUTINE anti_aliasing SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry + USE fields, ONLY: phi, moments_e, moments_i + IMPLICIT NONE IF ( contains_ky0 ) THEN ! Electron moments IF(KIN_E) THEN DO iz=izs,ize DO ij=ijgs_e,ijge_e DO ip=ipgs_e,ipge_e DO ikx=2,Nkx/2 !symmetry at ky = 0 moments_e( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_e( ip,ij,iky_0,Nkx+2-iky,iz, :)) END DO ! must be real at origin moments_e(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_e(ip,ij, iky_0,ikx_0,iz, :)) END DO END DO END DO ENDIF ! Ion moments DO iz=izs,ize DO ij=ijgs_i,ijge_i DO ip=ipgs_i,ipge_i DO ikx=2,Nkx/2 !symmetry at ky = 0 moments_i( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_i( ip,ij,iky_0,Nkx+2-ikx,iz, :)) END DO ! must be real at origin and top right moments_i(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_i(ip,ij, iky_0,ikx_0,iz, :)) END DO END DO END DO ! Phi DO ikx=2,Nkx/2 !symmetry at ky = 0 phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize) END DO ! must be real at origin phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize)) ENDIF END SUBROUTINE enforce_symmetry END SUBROUTINE stepon diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m index 3565c3d..5dc8efc 100644 --- a/wk/analysis_HeLaZ.m +++ b/wk/analysis_HeLaZ.m @@ -1,207 +1,207 @@ addpath(genpath([helazdir,'matlab'])) % ... add addpath(genpath([helazdir,'matlab/plot'])) % ... add addpath(genpath([helazdir,'matlab/compute'])) % ... add addpath(genpath([helazdir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [helazdir,'results/',outfile,'/']; MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax -JOBNUMMIN = 03; JOBNUMMAX = 20; +JOBNUMMIN = 00; JOBNUMMAX = 20; data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) options.TAVG_0 = 0.7*data.Ts3D(end); data.scale = (1/data.Nx/data.Ny)^2; options.TAVG_1 = data.Ts3D(end); % Averaging times duration options.NMVA = 1; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 1; fig = plot_radial_transport_and_spacetime(data,options); save_figure(data,fig,'.png') end if 0 %% statistical transport averaging options.T = [200 400]; fig = statistical_transport_averaging(data,options); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 1; % options.TIME = data.Ts5D(end-30:end); -options.TIME = data.Ts3D(1:end); +options.TIME = data.Ts3D(1:2:end); % options.TIME = [350:600]; data.EPS = 0.1; data.a = data.EPS * 2000; create_film(data,options,'.gif') end if 0 %% 2D snapshots % Options options.INTERP = 0; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'kxky'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 1; -options.TIME = [20 100 200 600 1900]; +options.TIME = [20 30 40 60 70]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 1; options.NAME = 'n_i'; options.PLANES = [1:10:80]; options.TIME = [200]; options.PLT_MTOPO = 0; data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) options.SPAR = linspace(-3,3,32)+(6/127/2); options.XPERP = linspace( 0,6,32); % options.SPAR = gene_data.vp'; % options.XPERP = gene_data.mu'; options.iz = 'avg'; options.T = [600]; options.PLT_FCT = 'contour'; options.ONED = 0; options.non_adiab = 0; options.SPECIE = 'i'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; options.NORMALIZED = 1; options.JOBNUM = 0; options.TIME = [300:500]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum options.TIME = 650:800; options.NORM =1; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'avg'; options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 'avg'; options.COMPY = 'avg'; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = 0:8000; options.NMA = 1; options.NMODES = 15; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos fig = plot_metric(data); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 15]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig,'.png') end diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index f87810f..f7de9f8 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,31 +1,30 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" helazdir = '/home/ahoffman/HeLaZ/'; % Directory of the simulation (from results) % if 1% Local results % outfile ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits % outfile ='shearless_cyclone/128x64x16x5x3_Dim_90'; % outfile ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; % outfile ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; % outfile ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS % outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; % outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; % outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; %% CBC % outfile ='shearless_cyclone/64x32x16x5x3_CBC_080'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; % outfile ='shearless_cyclone/64x32x16x5x3_CBC_120'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_080'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_100'; % outfile ='shearless_cyclone/64x32x16x9x5_CBC_120'; % outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% ZPINCH outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3'; - run analysis_HeLaZ