diff --git a/matlab/fort.90 b/matlab/fort.90
new file mode 100644
index 0000000000000000000000000000000000000000..321aad3bcbfa44db497c042526803dec63c6cc90
--- /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 50efa00a92690ac936f8667c480a49e5c358fcd9..c49fc47eff29195904ab50932b0a2e57a89af892 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -4,21 +4,26 @@ subroutine auxval
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 e71b4d03b2e5001927fb57b172ea5237dc367078..788a3aa816056285e354e64f4ddcf6227c5b7d49 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -20,9 +20,6 @@ MODULE basic
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)
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 03978ff3e10b9244179a6cd7afe7a202f5067330..94aedfbaa6caaaa7811f02f5d0043c9158052d6c 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -10,109 +10,112 @@ MODULE fourier
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 a3c0286bdb3cca69e0b2ff41c8cfe9882c9ca72d..da4bf65557b6cc2cebeceac4c07d2eac1c3e6d16 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -37,6 +37,7 @@ MODULE grid
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
@@ -88,86 +89,79 @@ CONTAINS
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
diff --git a/src/inital.F90 b/src/inital.F90
index 7f50dbb8096918d19ab8d7f0619780f55b26b628..7fa64c9988b9909fd5b4c85187ab23816ef53384 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -61,37 +61,7 @@ SUBROUTINE init_moments
! 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
@@ -106,13 +76,13 @@ SUBROUTINE init_moments
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
@@ -125,7 +95,6 @@ SUBROUTINE init_moments
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
@@ -133,7 +102,7 @@ SUBROUTINE init_moments
END DO
END DO
- ENDIF
+ ! ENDIF
END SUBROUTINE init_moments
!******************************************************************************!
@@ -177,9 +146,9 @@ SUBROUTINE load_cp
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)
@@ -191,8 +160,8 @@ SUBROUTINE load_cp
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)
diff --git a/src/ppexit.F90 b/src/ppexit.F90
index 9995c91e9b0ac9291e6dfba9012c249f7b1c5902..b6b04cd91571d07a1c53bc9da81e48b69c829154 100644
--- a/src/ppexit.F90
+++ b/src/ppexit.F90
@@ -2,13 +2,13 @@ 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)
diff --git a/src/stepon.F90 b/src/stepon.F90
index 7f16f1adaab6cf1c521bdb9c3b2d8ff0598dfdfd..986435d707709162307e13cbff549f29e3592770 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -17,11 +17,9 @@ SUBROUTINE stepon
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)
@@ -37,27 +35,22 @@ SUBROUTINE stepon
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
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 20d1a9152edbc68fdaa44e170a7f495d25c5a3b6..a47bcc7aca8d4f1e2c9966129846091ec3c3dbdf 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -89,13 +89,12 @@ 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
@@ -173,7 +172,7 @@ 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])
@@ -294,7 +293,7 @@ fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position', [100, 100, 1200
% 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;
diff --git a/wk/fort.90 b/wk/fort.90
index 412c03cac1350770f03669d8484a2b92377a29a5..79b41df03a3338fb1e5c64f71e732825acb790ac 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,7 +1,7 @@
&BASIC
nrun = 100000000
- dt = 0.01
- tmax = 1
+ dt = 0.005
+ tmax = 80
RESTART = .false.
/
&GRID
@@ -9,25 +9,25 @@
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