diff --git a/matlab/fort.90 b/matlab/fort.90 new file mode 100644 index 0000000..321aad3 --- /dev/null +++ b/matlab/fort.90 @@ -0,0 +1,6 @@ +&GRID + nr = 64 + Lr = 10 + nz = 64 + Lz = 10 +/ diff --git a/src/auxval.F90 b/src/auxval.F90 index 50efa00..c49fc47 100644 --- a/src/auxval.F90 +++ b/src/auxval.F90 @@ -1,24 +1,29 @@ subroutine auxval ! Set auxiliary values, at beginning of simulation USE basic USE grid USE array - USE fourier, ONLY: initialize_FFT + USE fourier, ONLY: init_grid_distr_and_plans use prec_const IMPLICIT NONE INTEGER :: irows,irowe, irow, icol IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ===' + CALL init_grid_distr_and_plans(Nr,Nz) + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + CALL set_krgrid + CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL set_kzgrid ! Distributed dimension + CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL set_pj + CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL memory ! Allocate memory for global arrays - - CALL initialize_FFT ! intialization of FFTW plans + CALL mpi_barrier(MPI_COMM_WORLD, ierr) END SUBROUTINE auxval diff --git a/src/basic_mod.F90 b/src/basic_mod.F90 index e71b4d0..788a3aa 100644 --- a/src/basic_mod.F90 +++ b/src/basic_mod.F90 @@ -1,281 +1,278 @@ MODULE basic ! Basic module for time dependent problems use, intrinsic :: iso_c_binding use prec_const IMPLICIT none ! INCLUDE 'fftw3-mpi.f03' INTEGER :: nrun = 1 ! Number of time steps to run real(dp) :: tmax = 100000.0 ! Maximum simulation time real(dp) :: dt = 1.0 ! Time step real(dp) :: time = 0 ! Current simulation time (Init from restart file) INTEGER :: jobnum = 0 ! Job number INTEGER :: step = 0 ! Calculation step of this run INTEGER :: cstep = 0 ! Current step number (Init from restart file) LOGICAL :: RESTART = .FALSE. ! Signal end of run LOGICAL :: nlend = .FALSE. ! Signal end of run INTEGER :: ierr ! flag for MPI error INTEGER :: my_id ! identification number of current process INTEGER :: num_procs ! number of MPI processes - integer(C_INTPTR_T) :: alloc_local ! total local data allocation size (mpi process) - integer(C_INTPTR_T) :: local_nkr ! local size of parallelized kz direction - integer(C_INTPTR_T) :: local_kr_start ! local start of parallelized kz direction INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose) INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose) INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose) INTEGER :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose) ! List of logical file units INTEGER :: lu_in = 90 ! File duplicated from STDIN INTEGER :: lu_job = 91 ! myjob file ! To measure computation time real :: start, finish real(dp) :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield, t0_step real(dp) :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield, t1_step real(dp) :: tc_rhs, tc_adv_field, tc_poisson, tc_Sapj, tc_diag, tc_checkfield, tc_step INTERFACE allocate_array MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4 MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5 MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4 MODULE PROCEDURE allocate_array_l1,allocate_array_l2,allocate_array_l3,allocate_array_l4 END INTERFACE allocate_array CONTAINS !================================================================================ SUBROUTINE basic_data ! Read basic data for input file use prec_const IMPLICIT NONE NAMELIST /BASIC/ nrun, dt, tmax, RESTART READ(lu_in,basic) !Init cumulative timers tc_rhs = 0.;tc_adv_field = 0.; tc_poisson = 0. tc_Sapj = 0.; tc_diag = 0.; tc_checkfield = 0. END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! Print date and time use prec_const IMPLICIT NONE CHARACTER(len=*) , INTENT(in) :: str CHARACTER(len=16) :: d, t, dat, time !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) time=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), time(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE display_h_min_s(time) real :: time integer :: days, hours, mins, secs days = FLOOR(time/24./3600.); hours= FLOOR(time/3600.); mins = FLOOR(time/60.); secs = FLOOR(time); IF ( days .GT. 0 ) THEN !display day h min s hours = (time/3600./24. - days) * 24 mins = (time/3600. - days*24. - hours) * 60 secs = (time/60. - days*24.*60 - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]' ELSEIF ( hours .GT. 0 ) THEN !display h min s mins = (time/3600. - hours) * 60 secs = (time/60. - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]' ELSEIF ( mins .GT. 0 ) THEN !display min s secs = (time/60. - mins) * 60 WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]' ELSE ! display s WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]' ENDIF END SUBROUTINE display_h_min_s !================================================================================ ! To allocate arrays of doubles, integers, etc. at run time SUBROUTINE allocate_array_dp1(a,is1,ie1) IMPLICIT NONE real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0.0_dp END SUBROUTINE allocate_array_dp1 SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2) IMPLICIT NONE real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=0.0_dp END SUBROUTINE allocate_array_dp2 SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0.0_dp END SUBROUTINE allocate_array_dp3 SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0.0_dp END SUBROUTINE allocate_array_dp4 SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=0.0_dp END SUBROUTINE allocate_array_dp5 !======================================== SUBROUTINE allocate_array_dc1(a,is1,ie1) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc1 SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc2 SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc3 SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc4 SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc5 !======================================== SUBROUTINE allocate_array_i1(a,is1,ie1) IMPLICIT NONE INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0 END SUBROUTINE allocate_array_i1 SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2) IMPLICIT NONE INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=0 END SUBROUTINE allocate_array_i2 SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0 END SUBROUTINE allocate_array_i3 SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0 END SUBROUTINE allocate_array_i4 SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=0 END SUBROUTINE allocate_array_i5 !======================================== SUBROUTINE allocate_array_l1(a,is1,ie1) IMPLICIT NONE LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=.false. END SUBROUTINE allocate_array_l1 SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2) IMPLICIT NONE LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=.false. END SUBROUTINE allocate_array_l2 SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=.false. END SUBROUTINE allocate_array_l3 SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=.false. END SUBROUTINE allocate_array_l4 SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=.false. END SUBROUTINE allocate_array_l5 END MODULE basic diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90 index 03978ff..94aedfb 100644 --- a/src/fourier_mod.F90 +++ b/src/fourier_mod.F90 @@ -1,118 +1,121 @@ MODULE fourier USE prec_const USE grid USE basic use, intrinsic :: iso_c_binding implicit none ! INCLUDE 'fftw3.f03' INCLUDE 'fftw3-mpi.f03' PRIVATE - type(C_PTR) :: planf, planb, real_ptr, cmpx_ptr - real(C_DOUBLE), pointer :: real_data_1(:,:), real_data_2(:,:) - complex(C_DOUBLE_complex), pointer :: cmpx_data(:,:) - integer (C_INTPTR_T) :: nx, ny, nh - PUBLIC :: initialize_FFT, finalize_FFT - PUBLIC :: convolve_2D_F2F + PUBLIC :: init_grid_distr_and_plans, convolve_2D_F2F, finalize_plans + real(C_DOUBLE), pointer :: real_data_f(:,:), real_data_g(:,:), real_data_c(:,:) + complex(C_DOUBLE_complex), pointer :: cmpx_data_f(:,:), cmpx_data_g(:,:), cmpx_data_c(:,:) + type(C_PTR) :: cdatar_f, cdatar_g, cdatar_c + type(C_PTR) :: cdatac_f, cdatac_g, cdatac_c + type(C_PTR) :: planf, planb + integer(C_INTPTR_T) :: i, ix, iy, alloc_local_1, alloc_local_2 + integer(C_INTPTR_T) :: NR_, NZ_ CONTAINS - SUBROUTINE initialize_FFT - IMPLICIT NONE - nx = Nr; ny = Nz; nh = ( nx / 2 ) + 1 + SUBROUTINE init_grid_distr_and_plans(Nr,Nz) + IMPLICIT NONE + + INTEGER, INTENT(IN) :: Nr,Nz + NR_ = Nr; NZ_ = Nz - IF ( my_id .EQ. 1)write(*,*) 'Initialize FFT..' - CALL fftw_mpi_init + ! Compute the room to allocate + alloc_local_1 = fftw_mpi_local_size_2d(NR_/2 + 1, NZ_, MPI_COMM_WORLD, local_nkr, local_nkr_offset) + ! Initalize pointers to this room + cdatac_f = fftw_alloc_complex(alloc_local_1) + cdatac_g = fftw_alloc_complex(alloc_local_1) + cdatac_c = fftw_alloc_complex(alloc_local_1) + ! Initalize the arrays with the rooms pointed + call c_f_pointer(cdatac_f, cmpx_data_f, [NZ_ ,local_nkr]) + call c_f_pointer(cdatac_g, cmpx_data_g, [NZ_ ,local_nkr]) + call c_f_pointer(cdatac_c, cmpx_data_c, [NZ_ ,local_nkr]) - IF ( my_id .EQ. 1)write(*,*) 'distribution of the array along kr..' - alloc_local = fftw_mpi_local_size_2d(nh, ny, MPI_COMM_WORLD, local_nkr, local_kr_start) - write(*,*) 'local_nkr', local_nkr, 'local_kr_start', local_kr_start + alloc_local_2 = fftw_mpi_local_size_2d(NZ_, NR_/2 + 1, MPI_COMM_WORLD, local_nz, local_nz_offset) - real_ptr = fftw_alloc_real(2 * alloc_local) - cmpx_ptr = fftw_alloc_complex(alloc_local) + cdatar_f = fftw_alloc_real(2*alloc_local_2) + cdatar_g = fftw_alloc_real(2*alloc_local_2) + cdatar_c = fftw_alloc_real(2*alloc_local_2) + call c_f_pointer(cdatar_f, real_data_f, [2*(NR_/2 + 1),local_nz]) + call c_f_pointer(cdatar_g, real_data_g, [2*(NR_/2 + 1),local_nz]) + call c_f_pointer(cdatar_c, real_data_c, [2*(NR_/2 + 1),local_nz]) - call c_f_pointer(real_ptr, real_data_1, [2*local_nkr, nx]) - call c_f_pointer(real_ptr, real_data_2, [2*local_nkr, nx]) - call c_f_pointer(cmpx_ptr, cmpx_data, [ local_nkr, nx]) + ! Plan Creation (out-of-place forward and backward FFT) + planf = fftw_mpi_plan_dft_r2c_2d(NZ_, NR_, real_data_f, cmpx_data_f, MPI_COMM_WORLD, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT)) + planb = fftw_mpi_plan_dft_c2r_2d(NZ_, NR_, cmpx_data_f, real_data_f, MPI_COMM_WORLD, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN)) - IF ( my_id .EQ. 1)write(*,*) 'setting FFT iFFT 2D plans..' - ! Backward FFT plan - planb = fftw_mpi_plan_dft_c2r_2d(ny, nx, cmpx_data, real_data_1, MPI_COMM_WORLD, & - FFTW_PATIENT) - ! Forward FFT plan - planf = fftw_mpi_plan_dft_r2c_2d(ny, nx, real_data_1, cmpx_data, MPI_COMM_WORLD, & - FFTW_PATIENT) - END SUBROUTINE initialize_FFT + if ((.not. c_associated(planf)) .OR. (.not. c_associated(planb))) then + write(*,*) "plan creation error!!" + stop + end if - SUBROUTINE finalize_FFT - IMPLICIT NONE + DO ix = 0,num_procs-1 + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. ix) print *, my_id,': alloc_local = ', alloc_local_1+alloc_local_2 + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! IF (my_id .EQ. 0) print *, my_id,': local_nz = ', local_nz, ' loc_nz_of. = ', local_nz_offset + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ENDDO - IF ( my_id .EQ. 0)write(*,*) 'finalize FFTW' - CALL dfftw_destroy_plan(planf) - CALL dfftw_destroy_plan(planb) - call fftw_mpi_cleanup - call fftw_free(real_ptr) - call fftw_free(cmpx_ptr) - END SUBROUTINE finalize_FFT + END SUBROUTINE init_grid_distr_and_plans !!! Convolution 2D Fourier to Fourier - ! - Compute the convolution using the convolution theorem + ! - Compute the convolution using the convolution theorem and MKL SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D ) IMPLICIT NONE - COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(IN) :: F_2D, G_2D ! input fields - COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(OUT) :: C_2D ! output convolutioned field - ! initialize data to some function my_function(i,j) - do ikr = ikrs,ikre - do ikz = ikzs,ikze - cmpx_data(ikr-local_kr_start, ikz) = F_2D(ikr, ikz) + COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(IN) :: F_2D, G_2D ! input fields + COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(OUT) :: C_2D ! output convolutioned field + + do ikr = ikrs, ikre + do ikz = ikzs, ikze + cmpx_data_f(ikz,ikr-local_nkr_offset) = F_2D(ikr,ikz) + cmpx_data_g(ikz,ikr-local_nkr_offset) = G_2D(ikr,ikz) end do end do - CALL mpi_barrier(MPI_COMM_WORLD, ierr) - ! 2D backward Fourier transform - ! write(*,*) my_id, 'iFFT(F) ..' - call fftw_mpi_execute_dft_c2r(planb, cmpx_data, real_data_1) + call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f) - ! initialize data to some function my_function(i,j) - do ikr = ikrs,ikre - do ikz = ikzs,ikze - cmpx_data(ikr-local_kr_start, ikz) = G_2D(ikr, ikz) - end do - end do - ! 2D inverse Fourier transform - ! write(*,*) my_id, 'iFFT(G) ..' - call fftw_mpi_execute_dft_c2r(planb, cmpx_data, real_data_2) - - ! Product in physical space - ! write(*,*) 'multply..' - real_data_1 = real_data_1 * real_data_2 - - ! 2D Fourier tranform - ! write(*,*) my_id, 'FFT(f*g) ..' - call fftw_mpi_execute_dft_r2c(planf, real_data_1, cmpx_data) - - ! COPY TO OUTPUT ARRAY - ! write(*,*) my_id, 'out ..' - do ikr = ikrs,ikre - do ikz = ikzs,ikze - cmpx_data(ikr-local_kr_start,ikz) = C_2D(ikr,ikz) + call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g) + + real_data_c = real_data_f/NZ_/NR_ * real_data_g/NZ_/NR_ + + call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c) + + ! Retrieve convolution in input format + do ikr = ikrs, ikre + do ikz = ikzs, ikze + C_2D(ikr,ikz) = cmpx_data_c(ikz,ikr-local_nkr_offset)*AA_r(ikr)*AA_z(ikz) end do end do - ! Anti aliasing (2/3 rule) - DO ikr = ikrs,ikre - DO ikz = ikzs,ikze - C_2D(ikr,ikz) = C_2D(ikr,ikz) * AA_r(ikr) * AA_z(ikz) - ENDDO - ENDDO END SUBROUTINE convolve_2D_F2F +SUBROUTINE finalize_plans + IMPLICIT NONE + + IF (my_id .EQ. 0) write(*,*) '..plan Destruction.' + call fftw_destroy_plan(planb) + call fftw_destroy_plan(planf) + call fftw_mpi_cleanup() + call fftw_free(cdatar_f) + call fftw_free(cdatar_g) + call fftw_free(cdatar_c) + call fftw_free(cdatac_f) + call fftw_free(cdatac_g) + call fftw_free(cdatac_c) +END SUBROUTINE finalize_plans + END MODULE fourier diff --git a/src/grid_mod.F90 b/src/grid_mod.F90 index a3c0286..da4bf65 100644 --- a/src/grid_mod.F90 +++ b/src/grid_mod.F90 @@ -1,221 +1,215 @@ 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 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 :: ikr - Nkr = Nr/2+1 ! Defined only on positive kr since fields are real - ! ! Start and END indices of grid - ikrs = my_id * Nkr/num_procs + 1 - ikre = (my_id+1) * Nkr/num_procs + SUBROUTINE set_krgrid + USE prec_const + IMPLICIT NONE + INTEGER :: i_ - IF(my_id .EQ. num_procs) ikre = Nkr + 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 - WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre + 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 - ! Grid spacings - IF (Lr .GT. 0) THEN + IF (my_id .EQ. num_procs-1) ikre = Nkr + !WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre + ! Grid spacings deltakr = 2._dp*PI/Lr - ELSE - deltakr = 1.0 - ENDIF - - ! Discretized kr positions ordered as dk*(0 1 2) - ALLOCATE(krarray(ikrs:ikre)) - DO ikr = ikrs,ikre - krarray(ikr) = REAL(ikr-1,dp) * deltakr - if (krarray(ikr) .EQ. 0) THEN - ikr_0 = ikr - ENDIF - 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) .GT. -two_third_krmax) .AND. (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 - USE model, ONLY : kr0KH, ikz0KH - IMPLICIT NONE - Nkz = Nz; - ! Start and END indices of grid - ikzs = 1 - ikze = Nkz - WRITE(*,*) 'ID = ',my_id,' ikzs = ', ikzs, ' ikze = ', ikze - ! Grid spacings - deltakz = 2._dp*PI/Lz - - ! 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 (kzarray(ikz) .EQ. 0) ikz_0 = ikz - if (ikz .EQ. Nz/2+1) kzarray(ikz) = -kzarray(ikz) - END DO - - IF (my_id .EQ. 1) WRITE(*,*) 'ID = ',my_id,' kz = ',kzarray - ! 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 - - ! Put kz0RT to the nearest grid point on kz - ikz0KH = NINT(kr0KH/deltakr)+1 - IF ( (ikz0KH .GE. ikzs) .AND. (ikz0KH .LE. ikze) ) kr0KH = kzarray(ikz0KH) - - END SUBROUTINE set_kzgrid + ! 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 + deltakz = 2._dp*PI/Lz + + ! 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 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 7f50dbb..7fa64c9 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,281 +1,250 @@ !******************************************************************************! !!!!!! 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 ENDIF !!!!!! Set phi !!!!!! IF (my_id .EQ. 1) WRITE(*,*) 'Init phi' CALL poisson !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!! IF ( NON_LIN .OR. (A0KH .NE. 0)) 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) - - IF ( only_Na00 ) THEN ! Spike initialization on density only - - DO ikr=ikrs,ikre - DO ikz=ikzs,ikze - moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) - moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) - END DO - END DO - - DO ikz=2,Nkz/2 !symmetry at kr = 0 - CALL RANDOM_NUMBER(noise) - moments_e( 1,1,1,ikz, :) = moments_e( 1,1,1,Nkz+2-ikz, :) - moments_i( 1,1,1,ikz, :) = moments_i( 1,1,1,Nkz+2-ikz, :) - END DO - - ELSE - !**** Gaussian initialization (Hakim 2017) ********************************* - ! sigma = 5._dp ! Gaussian sigma - ! gain = 0.5_dp ! Gaussian mean - ! kz_shift = 0.5_dp ! Gaussian z shift - ! moments_i = 0; moments_e = 0; - ! DO ikr=ikrs,ikre - ! kr = krarray(ikr) - ! DO ikz=ikzs,ikze - ! kz = kzarray(ikz) - ! moments_i( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp) - ! moments_e( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp) - ! END DO - ! END DO + 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 - CALL RANDOM_NUMBER(noise) 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 - CALL RANDOM_NUMBER(noise) moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :) END DO ENDIF END DO END DO - ENDIF + ! 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' WRITE(*,'(3x,a)') "Resume from previous run" CALL openf(rstfile, fidrst) IF (isgroup(fidrst,'/Basic/moments_e')) THEN 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_ 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 ! Read state of system from restart file - CALL getarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0) + 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/moments_i", n_ - CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0) + CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),pardim=3) ELSE CALL getatt(fidrst, '/Basic', 'cstep', cstep) CALL getatt(fidrst, '/Basic', 'time', time) CALL getatt(fidrst, '/Basic', 'jobnum', jobnum) jobnum = jobnum+1 CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d) CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d) iframe2d = iframe2d-1; iframe5d = iframe5d-1 ! Read state of system from restart file - CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0) - CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0) + CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),pardim=3) + CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),pardim=3) ENDIF CALL closef(fidrst) 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'); 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'); 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/ppexit.F90 b/src/ppexit.F90 index 9995c91..b6b04cd 100644 --- a/src/ppexit.F90 +++ b/src/ppexit.F90 @@ -1,15 +1,15 @@ SUBROUTINE ppexit ! Exit parallel environment USE basic - USE fourier, ONLY : finalize_FFT + USE fourier, ONLY : finalize_plans use prec_const IMPLICIT NONE - CALL finalize_FFT + CALL finalize_plans CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL mpi_finalize(ierr) END SUBROUTINE ppexit diff --git a/src/stepon.F90 b/src/stepon.F90 index 7f16f1a..986435d 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,110 +1,103 @@ 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 - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme) CALL advance_time_level - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! 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 - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Execution time end CALL cpu_time(t1_adv_field) tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field) - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Update electrostatic potential CALL poisson - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Update nonlinear term IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN CALL compute_Sapj ENDIF - CALL mpi_barrier(MPI_COMM_WORLD, ierr) !(The two routines above are called in inital for t=0) CALL checkfield_all() CALL mpi_barrier(MPI_COMM_WORLD, ierr) - 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 diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index 20d1a91..a47bcc7 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,379 +1,378 @@ %% Load results % JOBNUM = 0; load_results; % JOBNUM = 1; load_results; compile_results %% Retrieving max polynomial degree and sampling info Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe); Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi); Ns5D = numel(Ts5D); Ns2D = numel(Ts2D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts2D = Ts2D'; if strcmp(OUTPUTS.write_non_lin,'.true.') Si00 = squeeze(Sipj(1,1,:,:,:)); Se00 = squeeze(Sepj(1,1,:,:,:)); end %% Build grids Nkr = numel(kr); Nkz = numel(kz); [KZ,KR] = meshgrid(kz,kr); Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz); dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1); Lk = max(Lkr,Lkz); dr = 2*pi/Lk; dz = 2*pi/Lk; Nr = max(Nkr,Nkz); Nz = Nr; r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r); z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z); [ZZ,RR] = meshgrid(z,r); %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nr,Nz,Ns2D); ni00 = zeros(Nr,Nz,Ns2D); si00 = zeros(Nr,Nz,Ns5D); phi = zeros(Nr,Nz,Ns2D); drphi = zeros(Nr,Nz,Ns2D); dzphi = zeros(Nr,Nz,Ns2D); for it = 1:numel(Ts2D) NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz))); ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz))); phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz))); drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz))); dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz))); end if strcmp(OUTPUTS.write_non_lin,'.true.') for it = 1:numel(Ts5D) si00(:,:,it) = real(fftshift(ifft2(squeeze(Si00(:,:,it)),Nr,Nz))); end end % Post processing disp('- post processing') E_pot = zeros(1,Ns2D); % Potential energy n^2 E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2 ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi| Flux_ri = zeros(1,Ns2D); % Particle flux Gamma = Flux_zi = zeros(1,Ns2D); % Particle flux Gamma = Flux_re = zeros(1,Ns2D); % Particle flux Gamma = Flux_ze = zeros(1,Ns2D); % Particle flux Gamma = Ne_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj Ni_norm = zeros(Npi,Nji,Ns5D);% . if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj Si_norm = zeros(Npi,Nji,Ns5D);% . Sne00_norm = zeros(1,Ns2D); % Time evol. of the amp of e nonlin term Sni00_norm = zeros(1,Ns2D); % end Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2; for it = 1:numel(Ts2D) % Loop over 2D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); E_pot(it) = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id E_kin(it) = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr; ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz)))); Flux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; Flux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; end E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); dEdt = diff(E_pot+E_kin)./dt2D; for it = 1:numel(Ts5D) % Loop over 5D arrays - NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz; if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz; - Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz; Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz; + Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz; Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz; end end %% Compute growth rate disp('- growth rate') tend = Ts2D(end); tstart = 0.6*tend; g_ = zeros(Nkr,Nkz); [~,ikr0KH] = min(abs(kr-KR0KH)); for ikr = 1:Nkr for ikz = 1:Nkz g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend); end end % gkr0kz_Ni00 = max(real(g_(:,:)),[],1); gkr0kz_Ni00 = real(g_(ikr0KH,:)); %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 0 %% Time evolutions fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)]; set(gcf, 'Position', [100, 100, 900, 800]) subplot(221); for ip = 1:Npe for ij = 1:Nje plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$']; semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); subplot(222) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$']; semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); subplot(223) plot(Ts2D,Flux_ri,'-','DisplayName','$\Gamma_{ri}$'); hold on; plot(Ts2D,Flux_zi,'-','DisplayName','$\Gamma_{zi}$'); hold on; plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$') plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$') grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma$'); %legend('show'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$S_i^{',num2str(ip-1),num2str(ij-1),'}$']; semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on; end end grid on; xlabel('$t c_s/R$'); ylabel('$\sum_{k_r,k_z}|S_i^{pj}|$'); %legend('show'); else %% Growth rate subplot(224) [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz(1:Nz/2),gkr0kz_Ni00(1:Nz/2),... 'DisplayName',['J = ',num2str(JMAXI)]); xlabel('$k_z$'); ylabel('$\gamma_{Ni}$'); xlim([0. 1.0]); %ylim([0. 0.04]); end suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]); save_figure end %% if 0 %% Photomaton : real space -tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf)); +tf = 0; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['photo_real',sprintf('_t=%.0f',Ts2D(it))]; set(gcf, 'Position', [100, 100, 1500, 500]) subplot(131); plt = @(x) (((x))); pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$'); legend('$n_i$'); subplot(132); plt = @(x) ((x)); DATA = plt(ni00(:,:,it))-plt(ne00(:,:,it)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); legend('$n_i-n_e$'); set(gca,'ytick',[]); subplot(133); plt = @(x) ((x)); DATA = plt(phi(:,:,it)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); set(gca,'ytick',[]); legend('$\phi$'); % if strcmp(OUTPUTS.write_non_lin,'.true.') % subplot(133); plt = @(x) fftshift((abs(x)),2); % pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) % xlabel('$r/\rho_s$'); legend('$|S_i^{00}|$'); set(gca,'ytick',[]) % end suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB), sprintf(', $t c_s/R=%.0f$',Ts2D(it))]); save_figure end if 1 %% Photomaton : real space % FIELD = ni00; FNAME = 'ni'; FIELD = phi; FNAME = 'phi'; tf = 60; [~,it1] = min(abs(Ts2D-tf)); tf = 65; [~,it2] = min(abs(Ts2D-tf)); tf = 200; [~,it3] = min(abs(Ts2D-tf)); tf = 400; [~,it4] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position', [100, 100, 1500, 400]) subplot(141); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it1)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it1))); subplot(142); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it2)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); subplot(143); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it3)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it3))); subplot(144); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it4)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it4))); % suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),... % ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']); save_figure end %% t0 = 0; skip_ = 1; DELAY = 0.01*skip_; FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D); %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% Density ion GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_i$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Density electron GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Density ion - electron GIFNAME = ['ni-ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ni00+ne00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_i-n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Phi GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1; FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% Density ion frequency GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% Density ion frequency @ kr = 0 GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.'; X = (kz); T = Ts2D; YMIN = -.1; YMAX = 1.1; XMIN = min(kz); XMAX = max(kz); FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r\rho_s$'; create_gif_1D end %% if 1 %% Space time diagramm (fig 11 Ivanov 2020) t0 = 1; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D)); fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position', [100, 100, 1200, 400]) subplot(211) plot(Ts2D(it0:it1),Flux_ri(it0:it1)); ylabel('$\Gamma_r$'); grid on % title(['$\eta=',num2str(ETAB),'\quad',... % '\nu_{',CONAME,'}=',num2str(NU),'$']) % legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) - ylim([-0.01,1.1*max(Flux_ri(it0:it1))]); + ylim([0,1.1*max(Flux_ri(it0:it1))]); subplot(212) [TY,TX] = meshgrid(r,Ts2D(it0:it1)); pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,it0:it1),2))'); set(pclr, 'edgecolor','none'); %colorbar; xlabel('$t c_s/R$'), ylabel('$r/\rho_s$') legend('$\langle\partial_r \phi\rangle_z$') save_figure end if 0 %% Flow skip = 10; plt = @(x) x(1:skip:end,1:skip:end); tf = 120; [~,it] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = 'space_time_drphi'; quiver(plt(RR),plt(ZZ),plt(dzphi(:,:,it)),plt(-drphi(:,:,it)),0); xlabel('$r$'), ylabel('$z$') title(sprintf('$v_E$, $t=%.0f c_s/R$',Ts2D(it))); end if 0 %% Growth rate fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)]; [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz/kr(ikr0KH),gkr0kz_Ni00/(kr(ikr0KH)*A0KH),'DisplayName',['J = ',num2str(JMAXI)]) xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$'); xlim([0. 1.0]); ylim([0. 0.6]); save_figure end %% if 0 %% Mode time evolution [~,ik00] = min(abs(kz)); [~,idk] = min(abs(kz-dkz)); [~,ik50] = min(abs(kz-0.5*KR0KH)); [~,ik75] = min(abs(kz-0.7*KR0KH)); [~,ik10] = min(abs(kz-1.00*KR0KH)); plt = @(x) abs(squeeze(x)); fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)]; semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik00)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(idk)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik50)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik75)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik10)/KR0KH)]); hold on xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show'); semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts2D-Ts2D(end))),'k--') plot(tstart*[1 1],ylim,'k-','LineWidth',0.5); plot(tend*[1 1],ylim,'k-','LineWidth',0.5); save_figure end %% if 0 %% Show frame in kspace tf = 20; [~,it] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['krkz_frame',sprintf('_%.2d',JOBNUM)]; subplot(221); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat\phi|$'); subplot(222); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_e^{00}|$'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$'); end save_figure end diff --git a/wk/fort.90 b/wk/fort.90 index 412c03c..79b41df 100644 --- a/wk/fort.90 +++ b/wk/fort.90 @@ -1,64 +1,64 @@ &BASIC nrun = 100000000 - dt = 0.01 - tmax = 1 + dt = 0.005 + tmax = 80 RESTART = .false. / &GRID pmaxe =2 jmaxe = 1 pmaxi = 2 jmaxi = 1 - Nr = 64 - Lr = 20 - Nz = 64 - Lz = 20 + Nr = 256 + Lr = 66 + Nz = 256 + Lz = 66 kpar = 0 CANCEL_ODD_P = .false. / &OUTPUT_PAR - nsave_0d = 1 + nsave_0d = 100 nsave_1d = 0 - nsave_2d = 1 - nsave_5d = 1 + nsave_2d = 100 + nsave_5d = 200 write_Na00 = .true. write_moments = .true. write_phi = .true. write_non_lin = .true. write_doubleprecision = .true. - resfile0 = '../results/test_parallel/64x32_L_20_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/outputs' - rstfile0 = '../results/test_parallel/64x32_L_20_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/checkpoint' + resfile0 = '../results/test_parallel/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/outputs' + rstfile0 = '../results/test_parallel/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/checkpoint' job2load = 0 / &MODEL_PAR ! Collisionality CO = -2 DK = .false. NON_LIN = .true. mu = 5e-06 nu = 0.01 tau_e = 1 tau_i = 1 sigma_e = 0.023338 sigma_i = 1 q_e = -1 q_i = 1 eta_n = 1 eta_T = 0 eta_B = 0.5 lambdaD = 0 kr0KH = 0 A0KH = 0 / &INITIAL_CON only_Na00 =.false. initback_moments =0 initnoise_moments =1e-05 iseed =42 selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5' eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' / &TIME_INTEGRATION_PAR numerical_scheme='RK4' / \ No newline at end of file