diff --git a/src/auxval.F90 b/src/auxval.F90 index 07f7acd..6c28104 100644 --- a/src/auxval.F90 +++ b/src/auxval.F90 @@ -1,32 +1,31 @@ subroutine auxval ! Set auxiliary values, at beginning of simulation USE basic USE grid USE array USE fourier, ONLY: init_grid_distr_and_plans, alloc_local_1, alloc_local_2 use prec_const IMPLICIT NONE INTEGER :: irows,irowe, irow, icol, i_ IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ===' CALL init_grid_distr_and_plans(Nr,Nz) - CALL set_krgrid + CALL set_krgrid ! MPI Distributed dimension - CALL set_kzgrid ! Distributed dimension + CALL set_kzgrid CALL set_pj CALL memory ! Allocate memory for global arrays DO i_ = 0,num_procs-1 CALL mpi_barrier(MPI_COMM_WORLD, ierr) IF (my_id .EQ. i_) WRITE (*,'(I2,A9,I3,A8,I3,A8,I3,A8,I3,A15,I6)') & i_,': ikrs = ', ikrs, ' ikre = ', ikre, 'ikzs = ', ikzs, ' ikze = ', ikze, & 'alloc_local = ', alloc_local_1+alloc_local_2 - CALL mpi_barrier(MPI_COMM_WORLD, ierr) ENDDO END SUBROUTINE auxval diff --git a/src/control.F90 b/src/control.F90 index d95e890..a3f251e 100644 --- a/src/control.F90 +++ b/src/control.F90 @@ -1,83 +1,88 @@ SUBROUTINE control ! Control the run use basic use prec_const IMPLICIT NONE CALL cpu_time(start) !________________________________________________________________________________ ! 1. Prologue ! 1.1 Initialize the parallel environment - IF (my_id .EQ. 1) WRITE(*,*) 'Initialize MPI...' + IF (my_id .EQ. 0) WRITE(*,*) 'Initialize MPI...' CALL ppinit - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...MPI initialized' + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...MPI initialized' CALL daytim('Start at ') ! 1.2 Define data specific to run - IF (my_id .EQ. 1) WRITE(*,*) 'Load basic data...' + IF (my_id .EQ. 0) WRITE(*,*) 'Load basic data...' CALL basic_data - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...basic data loaded.' + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...basic data loaded.' ! 1.3 Read input parameters from input file - IF (my_id .EQ. 1) WRITE(*,*) 'Read input parameters...' + IF (my_id .EQ. 0) WRITE(*,*) 'Read input parameters...' CALL readinputs - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...input parameters read' + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...input parameters read' ! 1.4 Set auxiliary values (allocate arrays, set grid, ...) - IF (my_id .EQ. 1) WRITE(*,*) 'Calculate auxval...' + IF (my_id .EQ. 0) WRITE(*,*) 'Calculate auxval...' CALL auxval - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...auxval calculated' + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...auxval calculated' ! 1.5 Initial conditions - IF (my_id .EQ. 1) WRITE(*,*) 'Create initial state...' + IF (my_id .EQ. 0) WRITE(*,*) 'Create initial state...' CALL inital - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...initial state created' + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial state created' ! 1.6 Initial diagnostics - IF (my_id .EQ. 1) WRITE(*,*) 'Initial diagnostics...' + IF (my_id .EQ. 0) WRITE(*,*) 'Initial diagnostics...' CALL diagnose(0) - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...initial diagnostics done' + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial diagnostics done' CALL FLUSH(stdout) + CALL mpi_barrier(MPI_COMM_WORLD, ierr) !________________________________________________________________________________ - IF (my_id .EQ. 1) WRITE(*,*) 'Time integration loop..' + IF (my_id .EQ. 0) WRITE(*,*) 'Time integration loop..' !________________________________________________________________________________ ! 2. Main loop DO CALL cpu_time(t0_step) ! Measuring time step = step + 1 cstep = cstep + 1 CALL stepon time = time + dt CALL tesend IF( nlend ) EXIT ! exit do loop CALL cpu_time(t0_diag) ! Measuring time CALL diagnose(step) CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) CALL cpu_time(t1_step); tc_step = tc_step + (t1_step - t0_step) END DO - CALL mpi_barrier(MPI_COMM_WORLD,ierr) - IF (my_id .EQ. 1) WRITE(*,'(a/)') '...time integration done' + IF (my_id .EQ. 0) WRITE(*,'(a/)') '...time integration done' !________________________________________________________________________________ ! 9. Epilogue CALL diagnose(-1) CALL endrun - IF (my_id .EQ. 1) CALL daytim('Done at ') + IF (my_id .EQ. 0) CALL daytim('Done at ') CALL ppexit END SUBROUTINE control diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90 index 528e21d..2634028 100644 --- a/src/fourier_mod.F90 +++ b/src/fourier_mod.F90 @@ -1,119 +1,124 @@ MODULE fourier USE prec_const USE grid USE basic use, intrinsic :: iso_c_binding implicit none ! INCLUDE 'fftw3.f03' INCLUDE 'fftw3-mpi.f03' PRIVATE PUBLIC :: init_grid_distr_and_plans, convolve_2D_F2F, finalize_plans real(C_DOUBLE), pointer, PUBLIC :: real_data_f(:,:), real_data_g(:,:), real_data_c(:,:) complex(C_DOUBLE_complex), pointer, PUBLIC :: 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) , PUBLIC :: planf, planb integer(C_INTPTR_T) :: i, ix, iy integer(C_INTPTR_T), PUBLIC :: alloc_local_1, alloc_local_2 - integer(C_INTPTR_T) :: NR_, NZ_ + integer(C_INTPTR_T) :: NR_, NZ_, NR_halved + + ! many plan data variables + integer(C_INTPTR_T) :: howmany=9 ! numer of eleemnt of the tensor + integer :: rank=3 ! rank of the transform + integer(C_INTPTR_T), dimension(2) :: fft_dims ! array containing data extent CONTAINS SUBROUTINE init_grid_distr_and_plans(Nr,Nz) IMPLICIT NONE INTEGER, INTENT(IN) :: Nr,Nz NR_ = Nr; NZ_ = Nz + NR_halved = NR_/2 + 1 + !! Complex arrays F, G ! 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) + alloc_local_1 = fftw_mpi_local_size_2d(NR_halved, NZ_, MPI_COMM_WORLD, local_nkr, local_nkr_offset) + ! alloc_local_1 = fftw_mpi_local_size_2d_many(2, NR_halved, 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]) - alloc_local_2 = fftw_mpi_local_size_2d(NZ_, NR_/2 + 1, MPI_COMM_WORLD, local_nz, local_nz_offset) - + !! Real arrays iFFT(F), iFFT(G) + ! Compute the room to allocate + alloc_local_2 = fftw_mpi_local_size_2d(NZ_, NR_halved, MPI_COMM_WORLD, local_nz, local_nz_offset) + ! alloc_local_2 = fftw_mpi_local_size_2d_many(2, NZ_, NR_halved, MPI_COMM_WORLD, local_nz, local_nz_offset) + ! Initalize pointers to this room 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) + ! Initalize the arrays with the rooms pointed 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]) ! 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 ((.not. c_associated(planf)) .OR. (.not. c_associated(planb))) then write(*,*) "plan creation error!!" stop end if - 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) - ENDDO - - END SUBROUTINE init_grid_distr_and_plans !!! Convolution 2D Fourier to Fourier ! - Compute the convolution using the convolution theorem and MKL SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D ) IMPLICIT NONE 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 fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f) 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 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/inital.F90 b/src/inital.F90 index 048769d..d3e3b8c 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' + IF (my_id .EQ. 0) WRITE(*,*) 'Init phi' CALL poisson ENDIF !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!! IF ( NON_LIN ) THEN; - IF (my_id .EQ. 1) WRITE(*,*) 'Init Sapj' + IF (my_id .EQ. 0) 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 !******************************************************************************!