Skip to content
Snippets Groups Projects
Commit b37744a1 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

compiles and seems to work fine. Tests ahead

parent 00337f61
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/fort.90 0 → 100644
+ 6
− 0
View file @ b37744a1
&GRID
nr = 64
Lr = 10
nz = 64
Lz = 10
/
src/auxval.F90
+ 8
− 3
View file @ b37744a1
......@@ -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
src/basic_mod.F90
+ 0
− 3
View file @ b37744a1
......@@ -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)
......
src/fourier_mod.F90
+ 79
− 76
View file @ b37744a1
......@@ -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
src/grid_mod.F90
+ 68
− 74
View file @ b37744a1
......@@ -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
......
src/inital.F90
+ 7
− 38
View file @ b37744a1
......@@ -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)
......
src/ppexit.F90
+ 2
− 2
View file @ b37744a1
......@@ -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)
......
src/stepon.F90
+ 0
− 7
View file @ b37744a1
......@@ -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
......
wk/analysis_2D.m
+ 3
− 4
View file @ b37744a1
......@@ -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;
......
wk/fort.90
+ 11
− 11
View file @ b37744a1
&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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment