diff --git a/src/grid_mod.F90 b/src/grid_mod.F90 index a6660db..77b6119 100644 --- a/src/grid_mod.F90 +++ b/src/grid_mod.F90 @@ -1,222 +1,216 @@ MODULE grid ! Grid module for spatial discretization USE prec_const USE basic IMPLICIT NONE PRIVATE ! GRID Namelist INTEGER, PUBLIC, PROTECTED :: pmaxe = 1 ! The maximal electron Hermite-moment computed INTEGER, PUBLIC, PROTECTED :: jmaxe = 1 ! The maximal electron Laguerre-moment computed INTEGER, PUBLIC, PROTECTED :: pmaxi = 1 ! The maximal ion Hermite-moment computed INTEGER, PUBLIC, PROTECTED :: jmaxi = 1 ! The maximal ion Laguerre-moment computed INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal ion Laguerre-moment computed INTEGER, PUBLIC, PROTECTED :: Nr = 16 ! Number of total internal grid points in r REAL(dp), PUBLIC, PROTECTED :: Lr = 1._dp ! horizontal length of the spatial box INTEGER, PUBLIC, PROTECTED :: Nz = 16 ! Number of total internal grid points in z REAL(dp), PUBLIC, PROTECTED :: Lz = 1._dp ! vertical length of the spatial box INTEGER, PUBLIC, PROTECTED :: Nkr = 8 ! Number of total internal grid points in kr REAL(dp), PUBLIC, PROTECTED :: Lkr = 1._dp ! horizontal length of the fourier box INTEGER, PUBLIC, PROTECTED :: Nkz = 16 ! Number of total internal grid points in kz REAL(dp), PUBLIC, PROTECTED :: Lkz = 1._dp ! vertical length of the fourier box REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component - LOGICAL, PUBLIC, PROTECTED :: CANCEL_ODD_P = .false. ! To cancel odd Hermite polynomials INTEGER, PUBLIC, PROTECTED :: pskip = 0 ! variable to skip p degrees or not ! For Orszag filter REAL(dp), PUBLIC, PROTECTED :: two_third_krmax REAL(dp), PUBLIC, PROTECTED :: two_third_kzmax ! 1D Antialiasing arrays (2/3 rule) REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_r REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_z ! Grids containing position in physical space REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: rarray REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray REAL(dp), PUBLIC, PROTECTED :: deltar, deltaz INTEGER, PUBLIC, PROTECTED :: irs, ire, izs, ize INTEGER, PUBLIC :: ir,iz ! counters integer(C_INTPTR_T), PUBLIC :: local_nkr, local_nkr_offset, local_nz, local_nz_offset ! Grids containing position in fourier space REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray REAL(dp), PUBLIC, PROTECTED :: deltakr, deltakz INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze, a INTEGER, PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin INTEGER, PUBLIC :: ikr, ikz, ip, ij ! counters ! Grid containing the polynomials degrees INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i ! Public Functions PUBLIC :: set_pj PUBLIC :: set_krgrid, set_kzgrid PUBLIC :: grid_readinputs, grid_outputinputs PUBLIC :: bare, bari CONTAINS SUBROUTINE set_pj USE prec_const IMPLICIT NONE INTEGER :: ip, ij - IF (CANCEL_ODD_P) THEN - pskip = 1 - ELSE - pskip = 0 - ENDIF ips_e = 1; ipe_e = pmaxe/(1+pskip) + 1 ips_i = 1; ipe_i = pmaxi/(1+pskip) + 1 ALLOCATE(parray_e(ips_e:ipe_e)) ALLOCATE(parray_i(ips_i:ipe_i)) - DO ip = ips_e,ipe_e; parray_e(ip) = (1+pskip)*(ip-1); END DO - DO ip = ips_i,ipe_i; parray_i(ip) = (1+pskip)*(ip-1); END DO + DO ip = ips_e,ipe_e; parray_e(ip) = (ip-1); END DO + DO ip = ips_i,ipe_i; parray_i(ip) = (ip-1); END DO ijs_e = 1; ije_e = jmaxe + 1 ijs_i = 1; ije_i = jmaxi + 1 ALLOCATE(jarray_e(ijs_e:ije_e)) ALLOCATE(jarray_i(ijs_i:ije_i)) DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO maxj = MAX(jmaxi, jmaxe) END SUBROUTINE set_pj SUBROUTINE set_krgrid USE prec_const IMPLICIT NONE INTEGER :: i_ Nkr = Nr/2+1 ! Defined only on positive kr since fields are real ! Start and END indices of grid ikrs = local_nkr_offset + 1 ikre = ikrs + local_nkr - 1 DO i_ = 0,num_procs-1 CALL mpi_barrier(MPI_COMM_WORLD, ierr) IF (my_id .EQ. i_) print *, i_,': ikrs = ', ikrs, ' ikre = ', ikre CALL mpi_barrier(MPI_COMM_WORLD, ierr) ENDDO IF (my_id .EQ. num_procs-1) ikre = Nkr !WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre ! Grid spacings IF (Lr .EQ. 0) THEN deltakr = 1._dp ELSE deltakr = 2._dp*PI/Lr ENDIF ! Discretized kr positions ordered as dk*(0 1 2 3) ALLOCATE(krarray(ikrs:ikre)) DO ikr = ikrs,ikre krarray(ikr) = REAL(ikr-1,dp) * deltakr IF (krarray(ikr) .EQ. 0) ikr_0 = ikr END DO ! Orszag 2/3 filter two_third_krmax = 2._dp/3._dp*deltakr*Nkr ALLOCATE(AA_r(ikrs:ikre)) DO ikr = ikrs,ikre IF ( (krarray(ikr) .LT. two_third_krmax) ) THEN AA_r(ikr) = 1._dp; ELSE AA_r(ikr) = 0._dp; ENDIF END DO END SUBROUTINE set_krgrid SUBROUTINE set_kzgrid USE prec_const IMPLICIT NONE INTEGER :: i_ Nkz = Nz; ! Start and END indices of grid ikzs = 1 ikze = Nkz ! Grid spacings IF (Lz .EQ. 0) THEN deltakz = 1._dp ELSE deltakz = 2._dp*PI/Lz ENDIF ! Discretized kz positions ordered as dk*(0 1 2 3 -2 -1) ALLOCATE(kzarray(ikzs:ikze)) DO ikz = ikzs,ikze kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz))) if (ikz .EQ. Nz/2+1) kzarray(ikz) = -kzarray(ikz) IF (kzarray(ikz) .EQ. 0) ikz_0 = ikz END DO ! Orszag 2/3 filter two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2); ALLOCATE(AA_z(ikzs:ikze)) DO ikz = ikzs,ikze IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN AA_z(ikz) = 1._dp; ELSE AA_z(ikz) = 0._dp; ENDIF END DO END SUBROUTINE set_kzgrid SUBROUTINE grid_readinputs ! Read the input parameters USE prec_const IMPLICIT NONE INTEGER :: lu_in = 90 ! File duplicated from STDIN NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, & - Nr, Lr, Nz, Lz, kpar, CANCEL_ODD_P + Nr, Lr, Nz, Lz, kpar READ(lu_in,grid) END SUBROUTINE grid_readinputs SUBROUTINE grid_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "pmaxe", pmaxe) CALL attach(fidres, TRIM(str), "jmaxe", jmaxe) CALL attach(fidres, TRIM(str), "pmaxi", pmaxi) CALL attach(fidres, TRIM(str), "jmaxi", jmaxi) CALL attach(fidres, TRIM(str), "nr", nr) CALL attach(fidres, TRIM(str), "Lr", Lr) CALL attach(fidres, TRIM(str), "nz", nz) CALL attach(fidres, TRIM(str), "Lz", Lz) CALL attach(fidres, TRIM(str), "nkr", nkr) CALL attach(fidres, TRIM(str), "Lkr", Lkr) CALL attach(fidres, TRIM(str), "nkz", nkz) CALL attach(fidres, TRIM(str), "Lkz", Lkz) END SUBROUTINE grid_outputinputs FUNCTION bare(p_,j_) IMPLICIT NONE INTEGER :: bare, p_, j_ bare = (jmaxe+1)*p_ + j_ + 1 END FUNCTION FUNCTION bari(p_,j_) IMPLICIT NONE INTEGER :: bari, p_, j_ bari = (jmaxi+1)*p_ + j_ + 1 END FUNCTION END MODULE grid diff --git a/src/inital.F90 b/src/inital.F90 index a5ad1ed..048769d 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,242 +1,242 @@ !******************************************************************************! !!!!!! initialize the moments and load/build coeff tables !******************************************************************************! SUBROUTINE inital USE basic USE model, ONLY : CO, NON_LIN USE prec_const USE coeff USE time_integration USE array, ONLY : Sepj,Sipj implicit none real :: start_init, end_init, time_estimation CALL set_updatetlevel(1) !!!!!! Set the moments arrays Nepj, Nipj !!!!!! ! WRITE(*,*) 'Init moments' IF ( RESTART ) THEN CALL load_cp ELSE CALL init_moments !!!!!! Set phi !!!!!! IF (my_id .EQ. 1) WRITE(*,*) 'Init phi' CALL poisson ENDIF !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!! - IF ( NON_LIN .OR. (A0KH .NE. 0)) THEN; + IF ( NON_LIN ) THEN; IF (my_id .EQ. 1) WRITE(*,*) 'Init Sapj' CALL compute_Sapj ! WRITE(*,*) 'Building Dnjs table' CALL build_dnjs_table ENDIF !!!!!! Load the full coulomb collision operator coefficients !!!!!! IF (CO .EQ. -1) THEN IF (my_id .EQ. 0) WRITE(*,*) '=== Load Full Coulomb matrix ===' CALL load_FC_mat IF (my_id .EQ. 0) WRITE(*,*) '..done' ENDIF END SUBROUTINE inital !******************************************************************************! !******************************************************************************! !!!!!!! Initialize the moments randomly !******************************************************************************! SUBROUTINE init_moments USE basic USE grid USE fields USE prec_const USE utility, ONLY: checkfield USE initial_par IMPLICIT NONE REAL(dp) :: noise REAL(dp) :: kr, kz, sigma, gain, kz_shift INTEGER, DIMENSION(12) :: iseedarr ! Seed random number generator iseedarr(:)=iseed CALL RANDOM_SEED(PUT=iseedarr+my_id) !**** Broad noise initialization ******************************************* DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e DO ikr=ikrs,ikre DO ikz=ikzs,ikze CALL RANDOM_NUMBER(noise) moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) END DO END DO IF ( ikrs .EQ. 1 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :) END DO ENDIF END DO END DO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i DO ikr=ikrs,ikre DO ikz=ikzs,ikze CALL RANDOM_NUMBER(noise) moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) END DO END DO IF ( ikrs .EQ. 1 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :) END DO ENDIF END DO END DO ! ENDIF END SUBROUTINE init_moments !******************************************************************************! !******************************************************************************! !!!!!!! Load moments from a previous save !******************************************************************************! SUBROUTINE load_cp USE basic USE futils, ONLY: openf, closef, getarr, getatt, isgroup, isdataset USE grid USE fields USE diagnostics_par USE time_integration IMPLICIT NONE INTEGER :: rank, sz_, n_ INTEGER :: dims(1) = (/0/) CHARACTER(LEN=50) :: dset_name WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5' IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Resume from previous run" CALL openf(rstfile, fidrst,mpicomm=MPI_COMM_WORLD) n_ = 0 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ DO WHILE (isdataset(fidrst, dset_name)) n_ = n_ + 1 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ ENDDO n_ = n_ - 1 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_ ! Read state of system from restart file WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_ CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),pardim=3) WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ CALL getarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),pardim=3) WRITE(dset_name, "(A, '/', i6.6)") "/Basic/phi", n_ CALL getarr(fidrst, dset_name, phi(ikrs:ikre,ikzs:ikze),pardim=1) - + ! Read time dependent attributes CALL getatt(fidrst, dset_name, 'cstep', cstep) CALL getatt(fidrst, dset_name, 'time', time) CALL getatt(fidrst, dset_name, 'jobnum', jobnum) jobnum = jobnum+1 CALL getatt(fidrst, dset_name, 'iframe2d',iframe2d) CALL getatt(fidrst, dset_name, 'iframe5d',iframe5d) iframe2d = iframe2d-1; iframe5d = iframe5d-1 CALL closef(fidrst) IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!" END SUBROUTINE load_cp !******************************************************************************! !******************************************************************************! !!!!!!! Build the Laguerre-Laguerre coupling coefficient table !******************************************************************************! SUBROUTINE build_dnjs_table USE basic USE array, Only : dnjs USE grid, Only : jmaxe, jmaxi USE coeff IMPLICIT NONE INTEGER :: in, ij, is, J INTEGER :: n_, j_, s_ J = max(jmaxe,jmaxi) DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry n_ = in - 1 DO ij = in,J+1 j_ = ij - 1 DO is = ij,J+1 s_ = is - 1 dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0)) ! By symmetry dnjs(in,is,ij) = dnjs(in,ij,is) dnjs(ij,in,is) = dnjs(in,ij,is) dnjs(ij,is,in) = dnjs(in,ij,is) dnjs(is,ij,in) = dnjs(in,ij,is) dnjs(is,in,ij) = dnjs(in,ij,is) ENDDO ENDDO ENDDO END SUBROUTINE !******************************************************************************! !******************************************************************************! !!!!!!! Load the Full coulomb coefficient table from COSOlver results !******************************************************************************! SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6 USE basic USE grid USE initial_par USE array use model USE futils, ONLY : openf, getarr, closef IMPLICIT NONE INTEGER :: ip2,ij2 INTEGER :: fid1, fid2, fid3, fid4 !!!!!!!!!!!! Electron matrices !!!!!!!!!!!! ! get the self electron colision matrix CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD); CALL getarr(fid1, '/Caapj/Ceepj', Ceepj) ! get array (moli format) CALL closef(fid1) ! get the Test and Back field electron ion collision matrix CALL openf(eimat_file,fid2, 'r', 'D'); CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT) CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF) CALL closef(fid2) !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!! ! get the self electron colision matrix CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD); IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj CALL getarr(fid3, '/Caapj/Ceepj', Ciipj) ! get array (moli format) ELSE CALL getarr(fid3, '/Caapj/Ciipj', Ciipj) ! get array (moli format) ENDIF CALL closef(fid3) ! get the Test and Back field ion electron collision matrix CALL openf(iemat_file, fid4, 'r', 'D'); CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT) CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF) CALL closef(fid4) END SUBROUTINE load_FC_mat !******************************************************************************! diff --git a/src/memory.F90 b/src/memory.F90 index c9fa38e..58b8b39 100644 --- a/src/memory.F90 +++ b/src/memory.F90 @@ -1,41 +1,41 @@ SUBROUTINE memory ! Allocate arrays (done dynamically otherwise size is unknown) USE array USE basic USE fields USE grid USE time_integration - USE model, ONLY: CO, NON_LIN, A0KH + USE model, ONLY: CO, NON_LIN USE prec_const IMPLICIT NONE ! Moments and moments rhs CALL allocate_array( moments_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) ! Electrostatic potential CALL allocate_array(phi, ikrs,ikre, ikzs,ikze) ! Collision matrix IF (CO .EQ. -1) THEN CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1)) CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1)) CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1)) CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1)) CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1)) CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1)) ENDIF ! Non linear terms and dnjs table - IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN + IF ( NON_LIN ) THEN CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze ) CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze ) CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1) ENDIF END SUBROUTINE memory diff --git a/src/model_mod.F90 b/src/model_mod.F90 index d1f64b5..58604cf 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -1,74 +1,71 @@ MODULE model ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator LOGICAL, PUBLIC, PROTECTED :: DK = 0 ! Drift kinetic model LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp ! Mass REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: eta_n = 1._dp ! Density gradient REAL(dp), PUBLIC, PROTECTED :: eta_T = 1._dp ! Temperature gradient REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Magnetic gradient REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length - REAL(dp), PUBLIC :: kr0KH = 0._dp ! background wave frequ. for Rayleigh Taylor instability - INTEGER, PUBLIC :: ikz0KH = 0 ! '' index - REAL(dp), PUBLIC, PROTECTED :: A0KH = 0._dp ! background amplitude PUBLIC :: model_readinputs, model_outputinputs CONTAINS SUBROUTINE model_readinputs ! Read the input parameters USE basic, ONLY : lu_in USE prec_const IMPLICIT NONE NAMELIST /MODEL_PAR/ CO, DK, NON_LIN, mu, nu, tau_e, tau_i, sigma_e, sigma_i, & - q_e, q_i, eta_n, eta_T, eta_B, lambdaD, kr0KH, A0KH + q_e, q_i, eta_n, eta_T, eta_B, lambdaD READ(lu_in,model_par) !WRITE(*,model_par) ! Collision Frequency Normalization ... to match fluid limit nu = nu*0.532_dp END SUBROUTINE model_readinputs SUBROUTINE model_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "CO", CO) CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN) CALL attach(fidres, TRIM(str), "nu", nu) CALL attach(fidres, TRIM(str), "tau_e", tau_e) CALL attach(fidres, TRIM(str), "tau_i", tau_i) CALL attach(fidres, TRIM(str), "sigma_e", sigma_e) CALL attach(fidres, TRIM(str), "sigma_i", sigma_i) CALL attach(fidres, TRIM(str), "q_e", q_e) CALL attach(fidres, TRIM(str), "q_i", q_i) CALL attach(fidres, TRIM(str), "eta_n", eta_n) CALL attach(fidres, TRIM(str), "eta_T", eta_T) CALL attach(fidres, TRIM(str), "eta_B", eta_B) CALL attach(fidres, TRIM(str), "lambdaD", lambdaD) END SUBROUTINE model_outputinputs END MODULE model diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90 index cb5591f..8e25ed8 100644 --- a/src/moments_eq_rhs.F90 +++ b/src/moments_eq_rhs.F90 @@ -1,474 +1,474 @@ SUBROUTINE moments_eq_rhs USE basic USE time_integration USE array USE fields USE grid USE model USE prec_const USE utility, ONLY : is_nan IMPLICIT NONE INTEGER :: ip2, ij2, p_int, j_int, p2_int, j2_int ! loops indices and polynom. degrees REAL(dp) :: p_dp, j_dp REAL(dp) :: kr, kz, kperp2 REAL(dp) :: taue_qe_etaB, taui_qi_etaB REAL(dp) :: sqrtTaue_qe, sqrtTaui_qi, qe_sigmae_sqrtTaue, qi_sigmai_sqrtTaui REAL(dp) :: kernelj, kerneljp1, kerneljm1, be_2, bi_2 ! Kernel functions and variable REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables REAL(dp) :: xNapj, xNapp1j, xNapm1j, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop REAL(dp) :: xphij, xphijp1, xphijm1, xphijpar ! ESpot. factors depending on the pj loop REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs COMPLEX(dp) :: test_nan REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency ! Measuring execution time CALL cpu_time(t0_rhs) !Precompute species dependant factors IF( q_e .NE. 0._dp ) THEN taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term ELSE taue_qe_etaB = 0._dp sqrtTaue_qe = 0._dp ENDIF taui_qi_etaB = tau_i/q_i * eta_B ! factor of the magnetic moment coupling sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i) sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp nu_e = nu ! electron-ion collision frequency (where already multiplied by 0.532) nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ. nu_ee = nu_e/SQRT2 ! e-e coll. frequ. nu_ie = nu*sigma_e**2 ! i-e coll. frequ. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Electrons moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices p_int= parray_e(ip) ! Hermite polynom. degree p_dp = REAL(p_int,dp) ! REAL of the Hermite degree ! N_e^{p+1,j} coeff xNapp1j = sqrtTaue_qe * SQRT(p_dp + 1) ! N_e^{p-1,j} coeff xNapm1j = sqrtTaue_qe * SQRT(p_dp) ! N_e^{p+2,j} coeff xNapp2j = taue_qe_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) ! N_e^{p-2,j} coeff xNapm2j = taue_qe_etaB * SQRT(p_dp * (p_dp - 1._dp)) factj = 1.0 ! Start of the recursive factorial jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices j_int= jarray_e(ij) ! Laguerre polynom. degree j_dp = REAL(j_int,dp) ! REAL of degree ! N_e^{p,j+1} coeff xNapjp1 = -taue_qe_etaB * (j_dp + 1._dp) ! N_e^{p,j-1} coeff xNapjm1 = -taue_qe_etaB * j_dp ! N_e^{pj} coeff xNapj = taue_qe_etaB * 2._dp*(p_dp + j_dp + 1._dp) !! Collision operator pj terms xCapj = -nu_e*(p_dp + 2._dp*j_dp) !DK Lenard-Bernstein basis ! Dougherty part IF ( CO .EQ. -2) THEN IF ((p_int .EQ. 2) .AND. (j_int .EQ. 0)) THEN ! kronecker pj20 xCa20 = nu_e * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((p_int .EQ. 0) .AND. (j_int .EQ. 1)) THEN ! kronecker pj01 xCa20 = -nu_e * SQRT2 * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((p_int .EQ. 1) .AND. (j_int .EQ. 0)) THEN ! kronecker pj10 xCa20 = 0._dp xCa01 = 0._dp xCa10 = nu_e ELSE xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp ENDIF ENDIF !! Electrostatic potential pj terms IF (p_int .EQ. 0) THEN ! kronecker p0 xphij = (eta_n + 2.*j_dp*eta_T - 2._dp*eta_B*(j_dp+1._dp) ) xphijp1 = -(eta_T - eta_B)*(j_dp+1._dp) xphijm1 = -(eta_T - eta_B)* j_dp xphijpar= 0._dp ELSE IF (p_int .EQ. 1) THEN ! kronecker p1 xphijpar = qe_sigmae_sqrtTaue xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp ELSE IF (p_int .EQ. 2) THEN ! kronecker p2 xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp ELSE xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp ENDIF ! Recursive factorial IF (j_dp .GT. 0) THEN factj = factj * j_dp ELSE factj = 1._dp ENDIF ! Loop on kspace krloope : DO ikr = ikrs,ikre kzloope : DO ikz = ikzs,ikze kr = krarray(ikr) ! Poloidal wavevector kz = kzarray(ikz) ! Toroidal wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector be_2 = kperp2 * sigmae2_taue_o2 ! Kernel argument !! Compute moments and mixing terms ! term propto N_e^{p,j} TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel) ! term propto N_e^{p+1,j} and kparallel - IF ( (ip+1 .LE. pmaxe+1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check + IF ( (ip+1 .LE. pmaxe+1) ) THEN ! OoB check TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel) ELSE TNapp1j = 0._dp ENDIF ! term propto N_e^{p-1,j} and kparallel - IF ( (ip-1 .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check + IF ( (ip-1 .GE. 1) ) THEN ! OoB check TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel) ELSE TNapm1j = 0._dp ENDIF ! term propto N_e^{p+2,j} IF (ip+(2-pskip) .LE. pmaxe+1) THEN ! OoB check TNapp2j = xNapp2j * moments_e(ip+(2-pskip),ij,ikr,ikz,updatetlevel) ELSE TNapp2j = 0._dp ENDIF ! term propto N_e^{p-2,j} IF (ip-(2-pskip) .GE. 1) THEN ! OoB check TNapm2j = xNapm2j * moments_e(ip-(2-pskip),ij,ikr,ikz,updatetlevel) ELSE TNapm2j = 0._dp ENDIF ! xterm propto N_e^{p,j+1} IF (ij+1 .LE. jmaxe+1) THEN ! OoB check TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel) ELSE TNapjp1 = 0._dp ENDIF ! term propto N_e^{p,j-1} IF (ij-1 .GE. 1) THEN ! OoB check TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel) ELSE TNapjm1 = 0._dp ENDIF !! Collision IF (CO .EQ. -2) THEN ! Dougherty Collision terms IF ( (pmaxe .GE. 2-pskip) ) THEN ! OoB check TColl20 = xCa20 * moments_e(3-pskip,1,ikr,ikz,updatetlevel) ELSE TColl20 = 0._dp ENDIF IF ( (jmaxe .GE. 1) ) THEN ! OoB check TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel) ELSE TColl01 = 0._dp ENDIF - IF ( (pmaxe .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB + odd number for Hermite degrees check + IF ( (pmaxe .GE. 1) ) THEN ! OoB + odd number for Hermite degrees check TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel) ELSE TColl10 = 0._dp ENDIF ! Total collisional term TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)& + TColl20 + TColl01 + TColl10 ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix) TColl = 0._dp ! Initialization ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms p2_int = parray_e(ip2) jloopee: DO ij2 = 1,jmaxe+1 j2_int = jarray_e(ij2) TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) & *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int)) & +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int))) ENDDO jloopee ENDDO ploopee ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms p2_int = parray_i(ip2) jloopei: DO ij2 = 1,jmaxi+1 j2_int = jarray_i(ij2) TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) & *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int))) END DO jloopei ENDDO ploopei ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel) ENDIF !! Electrical potential term IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2 kernelj = be_2**j_int * exp(-be_2)/factj kerneljp1 = kernelj * be_2 /(j_dp + 1._dp) IF ( be_2 .NE. 0 ) THEN kerneljm1 = kernelj * j_dp / be_2 ELSE kerneljm1 = 0._dp ENDIF Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz) ELSE Tphi = 0._dp ENDIF ! Sum of all linear terms moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = & -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)& -imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) & - mu*kperp2**2 * moments_e(ip,ij,ikr,ikz,updatetlevel) & + TColl ! Adding non linearity - IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN + IF ( NON_LIN ) THEN moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = & moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz) ENDIF END DO kzloope END DO krloope END DO jloope END DO ploope !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int= parray_i(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree ! N_i^{p+1,j} coeff xNapp1j = sqrtTaui_qi * SQRT(p_dp + 1) ! N_i^{p-1,j} coeff xNapm1j = sqrtTaui_qi * SQRT(p_dp) ! x N_i^{p+2,j} coeff xNapp2j = taui_qi_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) ! x N_i^{p-2,j} coeff xNapm2j = taui_qi_etaB * SQRT(p_dp * (p_dp - 1._dp)) factj = 1._dp ! Start of the recursive factorial jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 j_int= jarray_i(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree ! x N_i^{p,j+1} coeff xNapjp1 = -taui_qi_etaB * (j_dp + 1._dp) ! x N_i^{p,j-1} coeff xNapjm1 = -taui_qi_etaB * j_dp ! x N_i^{pj} coeff xNapj = taui_qi_etaB * 2._dp*(p_dp + j_dp + 1._dp) !! Collision operator pj terms xCapj = -nu_i*(p_dp + 2._dp*j_dp) !DK Lenard-Bernstein basis ! Dougherty part IF ( CO .EQ. -2) THEN IF ((p_int .EQ. 2) .AND. (j_int .EQ. 0)) THEN ! kronecker pj20 xCa20 = nu_i * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((p_int .EQ. 0) .AND. (j_int .EQ. 1)) THEN ! kronecker pj01 xCa20 = -nu_i * SQRT2 * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((p_int .EQ. 1) .AND. (j_int .EQ. 0)) THEN ! kronecker pj10 xCa20 = 0._dp xCa01 = 0._dp xCa10 = nu_i ELSE xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp ENDIF ENDIF !! Electrostatic potential pj terms IF (p_int .EQ. 0) THEN ! krokecker p0 xphij = (eta_n + 2._dp*j_dp*eta_T - 2._dp*eta_B*(j_dp+1._dp)) xphijp1 = -(eta_T - eta_B)*(j_dp+1._dp) xphijm1 = -(eta_T - eta_B)* j_dp xphijpar = 0._dp ELSE IF (p_int .EQ. 1) THEN ! kronecker p1 xphijpar = qi_sigmai_sqrtTaui xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp ELSE IF (p_int .EQ. 2) THEN !krokecker p2 xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar = 0._dp ELSE xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar = 0._dp ENDIF ! Recursive factorial IF (j_dp .GT. 0) THEN factj = factj * j_dp ELSE factj = 1._dp ENDIF ! Loop on kspace krloopi : DO ikr = ikrs,ikre kzloopi : DO ikz = ikzs,ikze kr = krarray(ikr) ! Poloidal wavevector kz = kzarray(ikz) ! Toroidal wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector bi_2 = kperp2 * sigmai2_taui_o2 ! Kernel argument !! Compute moments and mixing terms ! term propto N_i^{p,j} TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel) ! term propto N_i^{p+1,j} - IF ( (ip+1 .LE. pmaxi+1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check + IF ( (ip+1 .LE. pmaxi+1) ) THEN ! OoB check TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel) ELSE TNapp1j = 0._dp ENDIF ! term propto N_i^{p-1,j} - IF ( (ip-1 .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check + IF ( (ip-1 .GE. 1) ) THEN ! OoB check TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel) ELSE TNapm1j = 0._dp ENDIF ! term propto N_i^{p+2,j} IF (ip+(2-pskip) .LE. pmaxi+1) THEN ! OoB check TNapp2j = xNapp2j * moments_i(ip+(2-pskip),ij,ikr,ikz,updatetlevel) ELSE TNapp2j = 0._dp ENDIF ! term propto N_i^{p-2,j} IF (ip-(2-pskip) .GE. 1) THEN ! OoB check TNapm2j = xNapm2j * moments_i(ip-(2-pskip),ij,ikr,ikz,updatetlevel) ELSE TNapm2j = 0._dp ENDIF ! xterm propto N_i^{p,j+1} IF (ij+1 .LE. jmaxi+1) THEN ! OoB check TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel) ELSE TNapjp1 = 0._dp ENDIF ! term propto N_i^{p,j-1} IF (ij-1 .GE. 1) THEN ! OoB check TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel) ELSE TNapjm1 = 0._dp ENDIF !! Collision IF (CO .EQ. -2) THEN ! Dougherty Collision terms IF ( (pmaxi .GE. 2-pskip) ) THEN ! OoB check TColl20 = xCa20 * moments_i(3-pskip,1,ikr,ikz,updatetlevel) ELSE TColl20 = 0._dp ENDIF IF ( (jmaxi .GE. 1) ) THEN ! OoB check TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel) ELSE TColl01 = 0._dp ENDIF - IF ( (pmaxi .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check + IF ( (pmaxi .GE. 1) ) THEN ! OoB check TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel) ELSE TColl10 = 0._dp ENDIF ! Total collisional term TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)& + TColl20 + TColl01 + TColl10 ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb for ions (COSOlver matrix) !!! TColl = 0._dp ! Initialization ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms p2_int = parray_i(ip2) jloopii: DO ij2 = 1,jmaxi+1 j2_int = jarray_i(ij2) TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) & *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int)) & +nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int))) ENDDO jloopii ENDDO ploopii ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms p2_int = parray_e(ip2) jloopie: DO ij2 = 1,jmaxe+1 j2_int = jarray_e(ij2) TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) & *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int))) ENDDO jloopie ENDDO ploopie ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel) ENDIF !! Electrical potential term IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2 kernelj = bi_2**j_int * exp(-bi_2)/factj kerneljp1 = kernelj * bi_2 /(j_dp + 1._dp) IF ( bi_2 .NE. 0 ) THEN kerneljm1 = kernelj * j_dp / bi_2 ELSE kerneljm1 = 0._dp ENDIF Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz) ELSE Tphi = 0._dp ENDIF ! Sum of linear terms moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = & -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)& -imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) & - mu*kperp2**2 * moments_i(ip,ij,ikr,ikz,updatetlevel) & + TColl ! Adding non linearity - IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN + IF ( NON_LIN ) THEN moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = & moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz) ENDIF END DO kzloopi END DO krloopi END DO jloopi END DO ploopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs diff --git a/src/stepon.F90 b/src/stepon.F90 index f79f0be..221f8e6 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,104 +1,104 @@ SUBROUTINE stepon ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme) USE basic USE time_integration USE fields, ONLY: moments_e, moments_i, phi USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj USE grid USE advance_field_routine, ONLY: advance_time_level, advance_field USE model USE utility, ONLY: checkfield use prec_const IMPLICIT NONE INTEGER :: num_step DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4 ! Compute right hand side of moments hierarchy equation CALL moments_eq_rhs ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme) CALL advance_time_level ! Update the moments with the hierarchy RHS (step by step) ! Execution time start CALL cpu_time(t0_adv_field) DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e CALL advance_field(moments_e(ip,ij,:,:,:),moments_rhs_e(ip,ij,:,:,:)) ENDDO ENDDO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i CALL advance_field(moments_i(ip,ij,:,:,:),moments_rhs_i(ip,ij,:,:,:)) ENDDO ENDDO ! Execution time end CALL cpu_time(t1_adv_field) tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field) ! Update electrostatic potential CALL poisson ! Update nonlinear term - IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN + IF ( NON_LIN ) THEN CALL compute_Sapj ENDIF !(The two routines above are called in inital for t=0) CALL checkfield_all() IF( nlend ) EXIT ! exit do loop END DO CONTAINS SUBROUTINE checkfield_all ! Check all the fields for inf or nan ! Execution time start CALL cpu_time(t0_checkfield) IF ( ikrs .EQ. 1 ) CALL enforce_symetry() ! Enforcing symmetry on kr = 0 IF(.NOT.nlend) THEN nlend=nlend .or. checkfield(phi,' phi') DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e nlend=nlend .or. checkfield(moments_e(ip,ij,:,:,updatetlevel),' moments_e') ENDDO ENDDO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i nlend=nlend .or. checkfield(moments_i(ip,ij,:,:,updatetlevel),' moments_i') ENDDO ENDDO ENDIF ! Execution time end CALL cpu_time(t1_checkfield) tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield) END SUBROUTINE checkfield_all SUBROUTINE enforce_symetry DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :)) END DO END DO END DO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :)) END DO END DO END DO DO ikz=2,Nkz/2 !symmetry at kr = 0 phi(1,ikz) = phi(1,Nkz+2-ikz) END DO END SUBROUTINE enforce_symetry END SUBROUTINE stepon