diff --git a/src/auxval.F90 b/src/auxval.F90 index 7d6fd96..db695e6 100644 --- a/src/auxval.F90 +++ b/src/auxval.F90 @@ -1,81 +1,93 @@ subroutine auxval ! Set auxiliary values, at beginning of simulation USE basic USE grid USE array USE model USE fourier, ONLY: init_grid_distr_and_plans, alloc_local_1, alloc_local_2 use prec_const USE numerics USE geometry IMPLICIT NONE INTEGER :: irows,irowe, irow, icol, i_ IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ===' IF (LINEARITY .NE. 'linear') THEN write(*,*) 'FFTW3 y-grid distribution' CALL init_grid_distr_and_plans(Nx,Ny) ELSE CALL init_1Dgrid_distr write(*,*) 'Manual y-grid distribution' ENDIF ! Init the grids CALL set_pgrid ! parallel kin (MPI distributed) CALL set_jgrid ! perp kin CALL set_kxgrid ! radial modes (MPI distributed by FFTW) CALL set_kygrid ! azymuthal modes CALL set_zgrid ! field aligned angle IF ((my_id .EQ. 0) .AND. SG) WRITE(*,*) '--2 staggered z grids--' CALL memory ! Allocate memory for global arrays CALL eval_magnetic_geometry ! precompute coeff for lin equation CALL compute_lin_coeff ! precompute coeff for lin equation and geometry CALL evaluate_kernels ! precompute the kernels CALL evaluate_poisson_op ! precompute the kernels IF ( LINEARITY .NE. 'linear' ) THEN; CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs ENDIF !! Display parallel settings DO i_ = 0,num_procs-1 CALL mpi_barrier(MPI_COMM_WORLD, ierr) IF (my_id .EQ. i_) THEN IF (my_id .EQ. 0) WRITE(*,*) '' IF (my_id .EQ. 0) WRITE(*,*) '--------- Parallel environement ----------' IF (my_id .EQ. 0) WRITE(*,'(A12,I3)') 'n_procs ', num_procs IF (my_id .EQ. 0) WRITE(*,'(A12,I3,A14,I3,A14,I3)') 'num_procs_p = ', num_procs_p, ', num_procs_ky = ', num_procs_ky, ', num_procs_z = ', num_procs_z IF (my_id .EQ. 0) WRITE(*,*) '' WRITE(*,'(A9,I3,A10,I3,A10,I3,A9,I3)')& 'my_id = ', my_id, ', rank_p = ', rank_p, ', rank_ky = ', rank_ky,', rank_z = ', rank_z WRITE(*,'(A22,I3,A11,I3)')& ' ips_e = ', ips_e, ', ipe_e = ', ipe_e WRITE(*,'(A22,I3,A11,I3)')& ' ijs_e = ', ijs_e, ', ije_e = ', ije_e WRITE(*,'(A22,I3,A11,I3)')& ' ips_i = ', ips_i, ', ipe_i = ', ipe_i WRITE(*,'(A22,I3,A11,I3)')& ' ijs_i = ', ijs_i, ', ije_i = ', ije_i WRITE(*,'(A22,I3,A11,I3)')& ' ikxs = ', ikxs , ', ikxe = ', ikxe WRITE(*,'(A22,I3,A11,I3)')& ' ikys = ', ikys , ', ikye = ', ikye WRITE(*,'(A22,I3,A11,I3)')& ' izs = ', izs , ', ize = ', ize IF (my_id .NE. num_procs-1) WRITE (*,*) '' IF (my_id .EQ. num_procs-1) WRITE(*,*) '------------------------------------------' ENDIF ENDDO CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN + IF(KIN_E) & + write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi + write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi + ENDIF + DO ip = ips_i,ipe_i + DO ij = ijs_i,ije_i + IF((parray_i(ip)+2*jarray_i(ij) .LE. dmaxi) .AND. (rank_ky + rank_z .EQ. 0))& + print*, '(',parray_i(ip),',',jarray_i(ij),')' + ENDDO + ENDDO + END SUBROUTINE auxval diff --git a/src/closure_mod.F90 b/src/closure_mod.F90 index b316e44..5d48206 100644 --- a/src/closure_mod.F90 +++ b/src/closure_mod.F90 @@ -1,150 +1,146 @@ module closure ! Contains the routines to define closures USE basic USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, nu, KIN_E USE grid USE array, ONLY: kernel_e, kernel_i USE fields, ONLY: moments_e, moments_i USE time_integration, ONLY: updatetlevel IMPLICIT NONE PUBLIC :: apply_closure_model CONTAINS ! Positive Oob indices are approximated with a model SUBROUTINE apply_closure_model IMPLICIT NONE CALL cpu_time(t0_clos) IF (CLOS .EQ. 0) THEN ! zero truncation, An+1=0 for n+1>nmax only CALL ghosts_upper_truncation ELSEIF (CLOS .EQ. 1) THEN - ! Truncation at highest fully represented kinetic moment + ! zero truncation, An+1=0 for n+1>nmax only + CALL ghosts_upper_truncation + ! Additional truncation at highest fully represented kinetic moment ! e.g. Dmax = 3 means - ! all Napj s.t. p+2j <= 3 + ! only Napj s.t. p+2j <= 3 are evolved ! -> (p,j) allowed are (0,0),(1,0),(0,1),(2,0),(1,1),(3,0) - ! =>> Dmax is Pmax, condition is p+2j<=Pmax - DO iz = izs,ize - DO ikx = ikxs,ikxe - DO iky = ikys,ikye - IF(KIN_E) THEN - DO ip = ipgs_e,ipge_e - DO ij = ijgs_e,ijge_e - IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) & - moments_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp - ENDDO - ENDDO - ENDIF - DO ip = ipgs_i,ipge_i - DO ij = ijgs_i,ijge_i - IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) & - moments_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp - ENDDO - ENDDO - ENDDO + ! =>> Dmax = min(Pmax,Jmax+1) + IF(KIN_E) THEN + DO ip = ipgs_e,ipge_e + DO ij = ijgs_e,ijge_e + IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) & + moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp + ENDDO + ENDDO + ENDIF + DO ip = ipgs_i,ipge_i + DO ij = ijgs_i,ijge_i + IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) & + moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDDO ENDDO ! + ghosts truncation CALL ghosts_upper_truncation ELSE ERROR STOP '! Closure scheme not found !' ENDIF CALL ghosts_lower_truncation CALL cpu_time(t1_clos) tc_clos = tc_clos + (t1_clos - t0_clos) END SUBROUTINE apply_closure_model ! Positive Oob indices are approximated with a model SUBROUTINE ghosts_upper_truncation IMPLICIT NONE ! zero truncation, An+1=0 for n+1>nmax ! Electrons IF(KIN_E) THEN ! applies only for the process that has largest j index IF(ije_e .EQ. Jmaxe+1) THEN DO ip = ipgs_e,ipge_e moments_e(ip,ije_e+1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDDO ENDIF ! applies only for the process that has largest p index IF(ipe_e .EQ. Pmaxe+1) THEN DO ij = ijgs_e,ijge_e moments_e(ipe_e+1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp IF(deltape .EQ. 1) THEN ! Must truncate the second stencil moments_e(ipe_e+2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDIF ENDDO ENDIF ENDIF ! Ions ! applies only for the process that has largest j index IF(ije_i .EQ. Jmaxi+1) THEN DO ip = ipgs_i,ipge_i moments_i(ip,ije_i+1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDDO ENDIF ! applies only for the process that has largest p index IF(ipe_i .EQ. Pmaxi+1) THEN DO ij = ijgs_i,ijge_i moments_i(ipe_i+1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp IF(deltape .EQ. 1) THEN ! Must truncate the second stencil moments_i(ipe_i+2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDIF ENDDO ENDIF END SUBROUTINE ghosts_upper_truncation ! Negative OoB indices are 0 SUBROUTINE ghosts_lower_truncation IMPLICIT NONE ! zero truncation, An=0 for n<0 ! Electrons IF(KIN_E) THEN ! applies only for the process that has lowest j index IF(ijs_e .EQ. 1) THEN DO ip = ipgs_e,ipge_e moments_e(ip,ijs_e-1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDDO ENDIF ! applies only for the process that has lowest p index IF(ips_e .EQ. 1) THEN DO ij = ijgs_e,ijge_e moments_e(ips_e-1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp IF(deltape .EQ. 1) THEN ! Must truncate the second stencil moments_e(ips_e-2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDIF ENDDO ENDIF ENDIF ! Ions IF(ijs_i .EQ. 1) THEN ! applies only for the process that has lowest j index DO ip = ipgs_i,ipge_i moments_i(ip,ijs_i-1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDDO ENDIF ! applies only for the process that has lowest p index IF(ips_i .EQ. 1) THEN DO ij = ijgs_i,ijge_i moments_i(ips_i-1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp IF(deltape .EQ. 1) THEN ! Must truncate the second stencil moments_i(ips_i-2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp ENDIF ENDDO ENDIF END SUBROUTINE ghosts_lower_truncation END module closure diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index 2fe9694..d2e316f 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -1,256 +1,256 @@ MODULE moments_eq_rhs IMPLICIT NONE PUBLIC :: compute_moments_eq_rhs CONTAINS SUBROUTINE compute_moments_eq_rhs USE model, only: KIN_E IMPLICIT NONE IF(KIN_E) CALL moments_eq_rhs_e CALL moments_eq_rhs_i END SUBROUTINE compute_moments_eq_rhs !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array USE fields USE grid USE model USE prec_const USE collision use geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2, dzlnB_o_J COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: i_ky,phikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloope : DO iz = izs,ize kxloope : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector kyloope : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value ! Kinetic loops jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1 j_int = jarray_e(ij) ploope : DO ip = ips_e, ipe_e ! Hermite loop p_int = parray_e(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxe)) THEN - !! Compute moments mixing terms - Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp - ! Perpendicular dynamic - ! term propto n_e^{p,j} - Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz) - ! term propto n_e^{p+2,j} - Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz) - ! term propto n_e^{p-2,j} - Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz) - ! term propto n_e^{p,j+1} - Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) - ! term propto n_e^{p,j-1} - Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) - ! Parallel dynamic - ! ddz derivative for Landau damping term - Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & - + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz) - ! Mirror terms - Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz) - Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz) - Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) - Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) - ! Trapping terms - Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) - Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz) - Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) + !! Compute moments mixing terms + Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp + ! Perpendicular dynamic + ! term propto n_e^{p,j} + Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz) + ! term propto n_e^{p+2,j} + Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz) + ! term propto n_e^{p-2,j} + Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz) + ! term propto n_e^{p,j+1} + Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) + ! term propto n_e^{p,j-1} + Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) + ! Parallel dynamic + ! ddz derivative for Landau damping term + Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & + + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz) + ! Mirror terms + Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz) + Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz) + Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) + Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) + ! Trapping terms + Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) + Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz) + Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) - Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 - !! Electrical potential term - IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 - Tphi = (xphij_i (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & - + xphijp1_i(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & - + xphijm1_i(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz - ELSE - Tphi = 0._dp - ENDIF + Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 + !! Electrical potential term + IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 + Tphi = (xphij_i (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + + xphijp1_i(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + + xphijm1_i(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz + ELSE + Tphi = 0._dp + ENDIF - !! Sum of all RHS terms - moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & - ! Perpendicular magnetic gradient/curvature effects - - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& - ! Parallel coupling (Landau Damping) - - Tpar*gradz_coeff(iz,eo) & - ! Mirror term (parallel magnetic gradient) - - gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) & - ! Drives (density + temperature gradients) - - i_ky * Tphi & - ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) - - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & - ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" - + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) & - ! Collision term - + TColl_e(ip,ij,iky,ikx,iz) & - ! Nonlinear term - - Sepj(ip,ij,iky,ikx,iz) + !! Sum of all RHS terms + moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & + ! Perpendicular magnetic gradient/curvature effects + - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& + ! Parallel coupling (Landau Damping) + - Tpar*gradz_coeff(iz,eo) & + ! Mirror term (parallel magnetic gradient) + - gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) & + ! Drives (density + temperature gradients) + - i_ky * Tphi & + ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) + - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & + ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" + + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) & + ! Collision term + + TColl_e(ip,ij,iky,ikx,iz) & + ! Nonlinear term + - Sepj(ip,ij,iky,ikx,iz) ELSE moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploope END DO jloope END DO kyloope END DO kxloope END DO zloope ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_e !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_i USE basic USE time_integration, ONLY: updatetlevel USE array USE fields USE grid USE model USE prec_const USE collision USE geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1 COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi COMPLEX(dp) :: i_ky, phikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloopi : DO iz = izs,ize kxloopi : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector kyloopi : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value ! Kinetic loops jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 j_int = jarray_i(ij) ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int = parray_i(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxi)) THEN - !! Compute moments mixing terms - Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp - ! Perpendicular dynamic - ! term propto n_i^{p,j} - Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz) - ! term propto n_i^{p+2,j} - Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz) - ! term propto n_i^{p-2,j} - Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz) - ! term propto n_i^{p,j+1} - Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz) - ! term propto n_i^{p,j-1} - Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz) - ! Tperp - Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 - ! Parallel dynamic - ! ddz derivative for Landau damping term - Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) & - + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz) - ! Mirror terms - Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz) - Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz) - Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) - Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) - ! Trapping terms - Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) - Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz) - Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) + !! Compute moments mixing terms + Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp + ! Perpendicular dynamic + ! term propto n_i^{p,j} + Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz) + ! term propto n_i^{p+2,j} + Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz) + ! term propto n_i^{p-2,j} + Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz) + ! term propto n_i^{p,j+1} + Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz) + ! term propto n_i^{p,j-1} + Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz) + ! Tperp + Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 + ! Parallel dynamic + ! ddz derivative for Landau damping term + Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) & + + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz) + ! Mirror terms + Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz) + Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz) + Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) + Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) + ! Trapping terms + Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) + Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz) + Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) - Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 + Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 - !! Electrical potential term - IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 - Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & - + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & - + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz - ELSE - Tphi = 0._dp - ENDIF + !! Electrical potential term + IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 + Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz + ELSE + Tphi = 0._dp + ENDIF - !! Sum of all RHS terms - moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & - ! Perpendicular magnetic gradient/curvature effects - - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp & - ! Parallel coupling (Landau damping) - - gradz_coeff(iz,eo) * Tpar & - ! Mirror term (parallel magnetic gradient) - - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & - ! Drives (density + temperature gradients) - - i_ky * Tphi & - ! Numerical hyperdiffusion (totally artificial, for stability purpose) - - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & - ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" - + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) & - ! Collision term - + TColl_i(ip,ij,iky,ikx,iz)& - ! Nonlinear term - - Sipj(ip,ij,iky,ikx,iz) + !! Sum of all RHS terms + moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & + ! Perpendicular magnetic gradient/curvature effects + - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp & + ! Parallel coupling (Landau damping) + - gradz_coeff(iz,eo) * Tpar & + ! Mirror term (parallel magnetic gradient) + - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & + ! Drives (density + temperature gradients) + - i_ky * Tphi & + ! Numerical hyperdiffusion (totally artificial, for stability purpose) + - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & + ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" + + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) & + ! Collision term + + TColl_i(ip,ij,iky,ikx,iz)& + ! Nonlinear term + - Sipj(ip,ij,iky,ikx,iz) ELSE moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploopi END DO jloopi END DO kyloopi END DO kxloopi END DO zloopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i END MODULE moments_eq_rhs diff --git a/src/processing_mod.F90 b/src/processing_mod.F90 index d4997a5..d98334c 100644 --- a/src/processing_mod.F90 +++ b/src/processing_mod.F90 @@ -1,490 +1,490 @@ MODULE processing ! contains the Hermite-Laguerre collision operators. Solved using COSOlver. USE basic USE prec_const USE grid USE geometry USE utility USE calculus implicit none REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri REAL(dp), PUBLIC, PROTECTED :: hflux_x PUBLIC :: compute_nadiab_moments_z_gradients_and_interp PUBLIC :: compute_density, compute_upar, compute_uperp PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments PUBLIC :: compute_radial_ion_transport, compute_radial_heatflux CONTAINS ! 1D diagnostic to compute the average radial particle transport _r SUBROUTINE compute_radial_ion_transport USE fields, ONLY : moments_i, phi USE array, ONLY : kernel_i USE geometry, ONLY : Jacobian USE time_integration, ONLY : updatetlevel IMPLICIT NONE COMPLEX(dp) :: pflux_local, gflux_local, integral REAL(dp) :: ky_, buffer(1:2) INTEGER :: i_, world_rank, world_size, root COMPLEX(dp), DIMENSION(izgs:izge) :: integrant pflux_local = 0._dp ! particle flux gflux_local = 0._dp ! gyrocenter flux integrant = 0._dp ! auxiliary variable for z integration IF(ips_i .EQ. 1) THEN !! Gyro center flux (drift kinetic) DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO iz = izgs,izge integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_i(ip0_i,1,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz)) ENDDO ! Integrate over z call simpson_rule_z(integrant,integral) ! add contribution gflux_local = gflux_local + integral*iInt_Jacobian ENDDO ENDDO !! Particle flux DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant = 0._dp ! auxiliary variable for z integration DO ij = ijs_i, ije_i DO iz = izgs,izge integrant(iz) = integrant(iz) + & Jacobian(iz,0)*imagu * ky_ * kernel_i(ij,iky,ikx,iz,0) & *moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz)) ENDDO ENDDO ! Integrate over z call simpson_rule_z(integrant,integral) ! add contribution pflux_local = pflux_local + integral*iInt_Jacobian ENDDO ENDDO ENDIF buffer(1) = REAL(gflux_local) buffer(2) = REAL(pflux_local) root = 0 !Gather manually among the rank_p=0 processes and perform the sum gflux_ri = 0 pflux_ri = 0 IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) gflux_ri = gflux_ri + buffer(1) pflux_ri = pflux_ri + buffer(2) ENDDO ENDIF ELSE gflux_ri = gflux_local pflux_ri = pflux_local ENDIF ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri END SUBROUTINE compute_radial_ion_transport ! 1D diagnostic to compute the average radial particle transport _r SUBROUTINE compute_radial_heatflux USE fields, ONLY : moments_i, moments_e, phi USE array, ONLY : kernel_e, kernel_i USE geometry, ONLY : Jacobian USE time_integration, ONLY : updatetlevel USE model, ONLY : qe_taue, qi_taui, KIN_E IMPLICIT NONE COMPLEX(dp) :: hflux_local, integral REAL(dp) :: ky_, buffer(1:2), j_dp INTEGER :: i_, world_rank, world_size, root COMPLEX(dp), DIMENSION(izgs:izge) :: integrant ! charge density q_a n_a hflux_local = 0._dp ! particle flux IF(ips_i .EQ. 1 .AND. ips_e .EQ. 1) THEN ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant = 0._dp DO ij = ijs_i, ije_i DO iz = izgs,izge integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))& *(twothird * ( 2._dp*j_dp * kernel_i(ij ,iky,ikx,iz,0) & - (j_dp+1) * kernel_i(ij+1,iky,ikx,iz,0) & - j_dp * kernel_i(ij-1,iky,ikx,iz,0))& * (moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel)+qi_taui*kernel_i(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) & + SQRT2*onethird * kernel_i(ij,iky,ikx,iz,0) * moments_i(ip2_i,ij,iky,ikx,iz,updatetlevel)) ENDDO ENDDO IF(KIN_E) THEN DO ij = ijs_e, ije_e DO iz = izgs,izge integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))& *(twothird * ( 2._dp*j_dp * kernel_e(ij ,iky,ikx,iz,0) & - (j_dp+1) * kernel_e(ij+1,iky,ikx,iz,0) & - j_dp * kernel_e(ij-1,iky,ikx,iz,0))& * (moments_e(ip0_e,ij,iky,ikx,iz,updatetlevel)+qe_taue*kernel_e(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) & + SQRT2*onethird * kernel_e(ij,iky,ikx,iz,0) * moments_e(ip2_e,ij,iky,ikx,iz,updatetlevel)) ENDDO ENDDO ENDIF ! Integrate over z call simpson_rule_z(integrant,integral) hflux_local = hflux_local + integral*iInt_Jacobian ENDDO ENDDO buffer(2) = REAL(hflux_local) root = 0 !Gather manually among the rank_p=0 processes and perform the sum hflux_x = 0 IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) hflux_x = hflux_x + buffer(2) ENDDO ENDIF ELSE hflux_x = hflux_local ENDIF ENDIF END SUBROUTINE compute_radial_heatflux SUBROUTINE compute_nadiab_moments_z_gradients_and_interp ! evaluate the non-adiabatique ion moments ! ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0 ! USE fields, ONLY : moments_i, moments_e, phi USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, & ddz_nepj, ddz2_Nepj, interp_nepj,& ddz_nipj, ddz2_Nipj, interp_nipj USE time_integration, ONLY : updatetlevel - USE model, ONLY : qe_taue, qi_taui, KIN_E + USE model, ONLY : qe_taue, qi_taui, KIN_E, CLOS USE calculus, ONLY : grad_z, grad_z2, interp_z IMPLICIT NONE - INTEGER :: eo, p_int + INTEGER :: eo, p_int, j_int CALL cpu_time(t0_process) ! Electron non adiab moments IF(KIN_E) THEN DO ip=ipgs_e,ipge_e IF(parray_e(ip) .EQ. 0) THEN DO ij=ijgs_e,ijge_e nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) & + qe_taue*kernel_e(ij,:,:,:,0)*phi(:,:,:) ENDDO ELSE DO ij=ijgs_e,ijge_e nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) ENDDO ENDIF ENDDO ENDIF ! Ions non adiab moments DO ip=ipgs_i,ipge_i IF(parray_i(ip) .EQ. 0) THEN DO ij=ijgs_i,ijge_i nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) & + qi_taui*kernel_i(ij,:,:,:,0)*phi(:,:,:) ENDDO ELSE DO ij=ijgs_i,ijge_i nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) ENDDO ENDIF ENDDO !! Ensure to kill the moments too high if closue option is set to 1 IF(CLOS .EQ. 1) THEN IF(KIN_E) THEN DO ip=ipgs_e,ipge_e p_int = parray_e(ip) DO ij=ijgs_e,ijge_e j_int = jarray_e(ij) IF(p_int+2*j_int .GT. dmaxe) & nadiab_moments_e(ip,ij,:,:,:) = 0._dp ENDDO ENDDO ENDIF DO ip=ipgs_i,ipge_i p_int = parray_i(ip) DO ij=ijgs_i,ijge_i j_int = jarray_i(ij) IF(p_int+2*j_int .GT. dmaxi) & nadiab_moments_i(ip,ij,:,:,:) = 0._dp ENDDO ENDDO ENDIF !------------- INTERP AND GRADIENTS ALONG Z ---------------------------------- IF (KIN_E) THEN DO ikx = ikxs,ikxe DO iky = ikys,ikye DO ij = ijgs_e,ijge_e DO ip = ipgs_e,ipge_e p_int = parray_e(ip) eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid ! Compute z derivatives CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize)) ! Parallel hyperdiffusion CALL grad_z2(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nepj(ip,ij,iky,ikx,izs:ize)) ! Compute even odd grids interpolation CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize)) ENDDO ENDDO ENDDO ENDDO ENDIF DO ikx = ikxs,ikxe DO iky = ikys,ikye DO ij = ijgs_i,ijge_i DO ip = ipgs_i,ipge_i p_int = parray_i(ip) eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid ! Compute z first derivative CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize)) ! Parallel hyperdiffusion CALL grad_z2(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nipj(ip,ij,iky,ikx,izs:ize)) ! Compute even odd grids interpolation CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize)) ENDDO ENDDO ENDDO ENDDO ! Execution time end CALL cpu_time(t1_process) tc_process = tc_process + (t1_process - t0_process) END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp !_____________________________________________________________________________! !!!!! FLUID MOMENTS COMPUTATIONS !!!!! ! Compute the 2D particle density for electron and ions (sum over Laguerre) SUBROUTINE compute_density USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: dens IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron density dens = 0._dp DO ij = ijs_e, ije_e dens = dens + kernel_e(ij,iky,ikx,iz,0) * nadiab_moments_e(ip0_e,ij,iky,ikx,iz) ENDDO dens_e(iky,ikx,iz) = dens ENDIF ! ion density dens = 0._dp DO ij = ijs_i, ije_i dens = dens + kernel_i(ij,iky,ikx,iz,0) * nadiab_moments_i(ip0_e,ij,iky,ikx,iz) ENDDO dens_i(iky,ikx,iz) = dens ENDDO ENDDO ENDDO ENDIF ! IF(KIN_E)& ! CALL manual_3D_bcast(dens_e(ikxs:ikxe,ikys:ikye,izs:ize)) ! CALL manual_3D_bcast(dens_i(ikxs:ikxe,ikys:ikye,izs:ize)) END SUBROUTINE compute_density ! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre) SUBROUTINE compute_uperp USE array, ONLY : uper_e, uper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: uperp IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron uperp = 0._dp DO ij = ijs_e, ije_e uperp = uperp + kernel_e(ij,iky,ikx,iz,0) *& 0.5_dp*(nadiab_moments_e(ip0_e,ij,iky,ikx,iz) - nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz)) ENDDO uper_e(iky,ikx,iz) = uperp ENDIF ! ion uperp = 0._dp DO ij = ijs_i, ije_i uperp = uperp + kernel_i(ij,iky,ikx,iz,0) *& 0.5_dp*(nadiab_moments_i(ip0_i,ij,iky,ikx,iz) - nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz)) ENDDO uper_i(iky,ikx,iz) = uperp ENDDO ENDDO ENDDO ENDIF ! IF(KIN_E)& ! CALL manual_3D_bcast(uper_e(ikxs:ikxe,ikys:ikye,izs:ize)) ! CALL manual_3D_bcast(uper_i(ikxs:ikxe,ikys:ikye,izs:ize)) END SUBROUTINE compute_uperp ! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre) SUBROUTINE compute_upar USE array, ONLY : upar_e, upar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: upar IF ( CONTAINS_ip1_e .AND. CONTAINS_ip1_i ) THEN DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron upar = 0._dp DO ij = ijs_e, ije_e upar = upar + kernel_e(ij,iky,ikx,iz,1)*nadiab_moments_e(ip1_e,ij,iky,ikx,iz) ENDDO upar_e(iky,ikx,iz) = upar ENDIF ! ion upar = 0._dp DO ij = ijs_i, ije_i upar = upar + kernel_i(ij,iky,ikx,iz,1)*nadiab_moments_i(ip1_i,ij,iky,ikx,iz) ENDDO upar_i(iky,ikx,iz) = upar ENDDO ENDDO ENDDO ELSE IF(KIN_E)& upar_e = 0 upar_i = 0 ENDIF ! IF(KIN_E)& ! CALL manual_3D_bcast(upar_e(ikxs:ikxe,ikys:ikye,izs:ize)) ! CALL manual_3D_bcast(upar_i(ikxs:ikxe,ikys:ikye,izs:ize)) END SUBROUTINE compute_upar ! Compute the 2D particle temperature for electron and ions (sum over Laguerre) SUBROUTINE compute_tperp USE array, ONLY : Tper_e, Tper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE REAL(dp) :: j_dp COMPLEX(dp) :: Tperp IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. & CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN ! Loop to compute T = 1/3*(Tpar + 2Tperp) DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe ! electron temperature IF(KIN_E) THEN Tperp = 0._dp DO ij = ijs_e, ije_e j_dp = REAL(ij-1,dp) Tperp = Tperp + kernel_e(ij,iky,ikx,iz,0)*& ((2_dp*j_dp+1)*nadiab_moments_e(ip0_e,ij ,iky,ikx,iz)& -j_dp *nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz)& -j_dp+1 *nadiab_moments_e(ip0_e,ij+1,iky,ikx,iz)) ENDDO Tper_e(iky,ikx,iz) = Tperp ENDIF ! ion temperature Tperp = 0._dp DO ij = ijs_i, ije_i j_dp = REAL(ij-1,dp) Tperp = Tperp + kernel_i(ij,iky,ikx,iz,0)*& ((2_dp*j_dp+1)*nadiab_moments_i(ip0_i,ij ,iky,ikx,iz)& -j_dp *nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz)& -j_dp+1 *nadiab_moments_i(ip0_i,ij+1,iky,ikx,iz)) ENDDO Tper_i(iky,ikx,iz) = Tperp ENDDO ENDDO ENDDO ENDIF ! IF(KIN_E)& ! CALL manual_3D_bcast(Tper_e(ikxs:ikxe,ikys:ikye,izs:ize)) ! CALL manual_3D_bcast(Tper_i(ikxs:ikxe,ikys:ikye,izs:ize)) END SUBROUTINE compute_Tperp ! Compute the 2D particle temperature for electron and ions (sum over Laguerre) SUBROUTINE compute_Tpar USE array, ONLY : Tpar_e, Tpar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE REAL(dp) :: j_dp COMPLEX(dp) :: tpar IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. & CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN ! Loop to compute T = 1/3*(Tpar + 2Tperp) DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe ! electron temperature IF(KIN_E) THEN Tpar = 0._dp DO ij = ijs_e, ije_e j_dp = REAL(ij-1,dp) Tpar = Tpar + kernel_e(ij,iky,ikx,iz,0)*& (SQRT2 * nadiab_moments_e(ip2_e,ij,iky,ikx,iz) & + nadiab_moments_e(ip0_e,ij,iky,ikx,iz)) ENDDO Tpar_e(iky,ikx,iz) = Tpar ENDIF ! ion temperature Tpar = 0._dp DO ij = ijs_i, ije_i j_dp = REAL(ij-1,dp) Tpar = Tpar + kernel_i(ij,iky,ikx,iz,0)*& (SQRT2 * nadiab_moments_i(ip2_i,ij,iky,ikx,iz) & + nadiab_moments_i(ip0_i,ij,iky,ikx,iz)) ENDDO Tpar_i(iky,ikx,iz) = Tpar ENDDO ENDDO ENDDO ENDIF ! IF(KIN_E)& ! CALL manual_3D_bcast(Tpar_e(ikxs:ikxe,ikys:ikye,izs:ize)) ! CALL manual_3D_bcast(Tpar_i(ikxs:ikxe,ikys:ikye,izs:ize)) END SUBROUTINE compute_Tpar ! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre) SUBROUTINE compute_fluid_moments USE array, ONLY : dens_e, Tpar_e, Tper_e, dens_i, Tpar_i, Tper_i, temp_e, temp_i USE model, ONLY : KIN_E IMPLICIT NONE CALL compute_density CALL compute_upar CALL compute_uperp CALL compute_Tpar CALL compute_Tperp ! Temperature temp_e = (Tpar_e - 2._dp * Tper_e)/3._dp - dens_e temp_i = (Tpar_i - 2._dp * Tper_i)/3._dp - dens_i END SUBROUTINE compute_fluid_moments END MODULE processing diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index 2add3a7..0f2cdef 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,194 +1,195 @@ 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]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 10; data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) options.TAVG_0 = 0.8*data.Ts3D(end); 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 = '\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) end if 0 %% statistical transport averaging options.T = [16000 17000]; fig = statistical_transport_averaging(data,options); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; -% options.NAME = '\phi'; +options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; -options.NAME = '\Gamma_x'; +% options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; -options.PLAN = 'kxky'; +options.PLAN = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; % options.TIME = dat.Ts5D; -options.TIME = 00:1:200; +options.TIME = 00:1:250; 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 = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 1; options.TIME = [150]; data.a = data.EPS * 1000; 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; options.TIME = [100 200]; options.PLT_MTOPO = 1; data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig) end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) options.SPAR = linspace(-3,3,64)+(6/127/2); options.XPERP = linspace( 0,6,64); % options.SPAR = vp'; % options.XPERP = mu'; options.Z = 1; -options.T = 5500; +options.T = 200; options.CTR = 1; options.ONED = 0; fig = plot_fa(data,options); save_figure(data,fig) end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; -options.P2J = 0; +options.P2J = 1; options.ST = 0; options.NORMALIZED = 0; fig = show_moments_spectrum(data,options); save_figure(data,fig) end if 0 %% Time averaged spectrum options.TIME = 1000:5000; options.NORM = 0; options.NAME = '\phi'; options.NAME = 'n_i'; options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 0; options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig) 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) end if 0 %% Mode evolution options.NORMALIZED = 1; options.K2PLOT = 1; options.TIME = 5:1:15; options.NMA = 1; options.NMODES = 5; options.iz = 1; fig = mode_growth_meter(data,options); save_figure(data,fig) end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 200; TAVG_1 = 3000; % 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) 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) end diff --git a/wk/analysis_header.m b/wk/analysis_header.m index b3cf14f..5188876 100644 --- a/wk/analysis_header.m +++ b/wk/analysis_header.m @@ -1,12 +1,15 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" helazdir = '/home/ahoffman/HeLaZ/'; % Directory of the simulation (from results) % if 1% Local results outfile =''; +outfile =''; +outfile =''; +outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0'; % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK'; % outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1'; % outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK'; -outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/'; +% outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/'; % outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/'; % outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1'; run analysis_3D