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

typos

parent 5228831e
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 638 additions and 635 deletions
src/calculus_mod.F90
+ 47
− 47
View file @ a0bebee6
...@@ -20,25 +20,25 @@ MODULE calculus ...@@ -20,25 +20,25 @@ MODULE calculus
CONTAINS CONTAINS
SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf) SUBROUTINE grad_z(target,local_nz,Ngz,inv_deltaz,f,ddzf)
implicit none implicit none
! Compute the periodic boundary condition 4 points centered finite differences ! Compute the periodic boundary condition 4 points centered finite differences
! formula among staggered grid or not. ! formula among staggered grid or not.
! not staggered : the derivative results must be on the same grid as the field ! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other ! staggered : the derivative is computed from a grid to the other
INTEGER, INTENT(IN) :: target, local_Nz, Ngz INTEGER, INTENT(IN) :: target, local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz REAL(dp), INTENT(IN) :: inv_deltaz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzf COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzf
INTEGER :: iz INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
SELECT CASE(TARGET) SELECT CASE(TARGET)
CASE(1) CASE(1)
CALL grad_z_o2e(local_Nz,Ngz,inv_deltaz,f,ddzf) CALL grad_z_o2e(local_nz,Ngz,inv_deltaz,f,ddzf)
CASE(2) CASE(2)
CALL grad_z_e2o(local_Nz,Ngz,inv_deltaz,f,ddzf) CALL grad_z_e2o(local_nz,Ngz,inv_deltaz,f,ddzf)
CASE DEFAULT ! No staggered grid -> usual centered finite differences CASE DEFAULT ! No staggered grid -> usual centered finite differences
DO iz = 1,local_Nz DO iz = 1,local_nz
ddzf(iz) = dz_usu(-2)*f(iz ) + dz_usu(-1)*f(iz+1) & ddzf(iz) = dz_usu(-2)*f(iz ) + dz_usu(-1)*f(iz+1) &
+dz_usu( 0)*f(iz+2) & +dz_usu( 0)*f(iz+2) &
+dz_usu( 1)*f(iz+3) + dz_usu( 2)*f(iz+4) +dz_usu( 1)*f(iz+3) + dz_usu( 2)*f(iz+4)
...@@ -49,30 +49,30 @@ SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf) ...@@ -49,30 +49,30 @@ SUBROUTINE grad_z(target,local_Nz,Ngz,inv_deltaz,f,ddzf)
ddzf = 0._dp ddzf = 0._dp
ENDIF ENDIF
CONTAINS CONTAINS
SUBROUTINE grad_z_o2e(local_Nz,Ngz,inv_deltaz,fo,ddzfe) ! Paruta 2018 eq (27) SUBROUTINE grad_z_o2e(local_nz,Ngz,inv_deltaz,fo,ddzfe) ! Paruta 2018 eq (27)
! gives the gradient of a field from the odd grid to the even one ! gives the gradient of a field from the odd grid to the even one
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz REAL(dp), INTENT(IN) :: inv_deltaz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fo COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fo
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzfe ! COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfe !
INTEGER :: iz INTEGER :: iz
DO iz = 1,local_Nz DO iz = 1,local_nz
ddzfe(iz) = dz_o2e(-2)*fo(iz ) + dz_o2e(-1)*fo(iz+1) & ddzfe(iz) = dz_o2e(-2)*fo(iz ) + dz_o2e(-1)*fo(iz+1) &
+dz_o2e( 0)*fo(iz+2) + dz_o2e( 1)*fo(iz+3) +dz_o2e( 0)*fo(iz+2) + dz_o2e( 1)*fo(iz+3)
ENDDO ENDDO
ddzfe(:) = ddzfe(:) * inv_deltaz ddzfe(:) = ddzfe(:) * inv_deltaz
END SUBROUTINE grad_z_o2e END SUBROUTINE grad_z_o2e
SUBROUTINE grad_z_e2o(local_Nz,Ngz,inv_deltaz,fe,ddzfo) ! n2v for Paruta 2018 eq (28) SUBROUTINE grad_z_e2o(local_nz,Ngz,inv_deltaz,fe,ddzfo) ! n2v for Paruta 2018 eq (28)
! gives the gradient of a field from the even grid to the odd one ! gives the gradient of a field from the even grid to the odd one
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz REAL(dp), INTENT(IN) :: inv_deltaz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fe COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fe
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddzfo COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddzfo
INTEGER :: iz INTEGER :: iz
DO iz = 1,local_Nz DO iz = 1,local_nz
ddzfo(iz) = dz_e2o(-1)*fe(iz+1) + dz_e2o(0)*fe(iz+2) & ddzfo(iz) = dz_e2o(-1)*fe(iz+1) + dz_e2o(0)*fe(iz+2) &
+dz_e2o( 1)*fe(iz+3) + dz_e2o(2)*fe(iz+4) +dz_e2o( 1)*fe(iz+3) + dz_e2o(2)*fe(iz+4)
ENDDO ENDDO
...@@ -80,16 +80,16 @@ CONTAINS ...@@ -80,16 +80,16 @@ CONTAINS
END SUBROUTINE grad_z_e2o END SUBROUTINE grad_z_e2o
END SUBROUTINE grad_z END SUBROUTINE grad_z
SUBROUTINE grad_z2(local_Nz,Ngz,inv_deltaz,f,ddz2f) SUBROUTINE grad_z2(local_nz,Ngz,inv_deltaz,f,ddz2f)
! Compute the second order fourth derivative for periodic boundary condition ! Compute the second order fourth derivative for periodic boundary condition
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz REAL(dp), INTENT(IN) :: inv_deltaz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddz2f COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz2f
INTEGER :: iz INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_Nz DO iz = 1,local_nz
ddz2f(iz) = dz2_usu(-2)*f(iz ) + dz2_usu(-1)*f(iz+1) & ddz2f(iz) = dz2_usu(-2)*f(iz ) + dz2_usu(-1)*f(iz+1) &
+dz2_usu( 0)*f(iz+2)& +dz2_usu( 0)*f(iz+2)&
+dz2_usu( 1)*f(iz+3) + dz2_usu( 2)*f(iz+4) +dz2_usu( 1)*f(iz+3) + dz2_usu( 2)*f(iz+4)
...@@ -101,16 +101,16 @@ SUBROUTINE grad_z2(local_Nz,Ngz,inv_deltaz,f,ddz2f) ...@@ -101,16 +101,16 @@ SUBROUTINE grad_z2(local_Nz,Ngz,inv_deltaz,f,ddz2f)
END SUBROUTINE grad_z2 END SUBROUTINE grad_z2
SUBROUTINE grad_z4(local_Nz,Ngz,inv_deltaz,f,ddz4f) SUBROUTINE grad_z4(local_nz,Ngz,inv_deltaz,f,ddz4f)
! Compute the second order fourth derivative for periodic boundary condition ! Compute the second order fourth derivative for periodic boundary condition
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
REAL(dp), INTENT(IN) :: inv_deltaz REAL(dp), INTENT(IN) :: inv_deltaz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: ddz4f COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: ddz4f
INTEGER :: iz INTEGER :: iz
IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid IF(Ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = 1,local_Nz DO iz = 1,local_nz
ddz4f(iz) = dz4_usu(-2)*f(iz ) + dz4_usu(-1)*f(iz+1) & ddz4f(iz) = dz4_usu(-2)*f(iz ) + dz4_usu(-1)*f(iz+1) &
+dz4_usu( 0)*f(iz+2)& +dz4_usu( 0)*f(iz+2)&
+dz4_usu( 1)*f(iz+3) + dz4_usu( 2)*f(iz+4) +dz4_usu( 1)*f(iz+3) + dz4_usu( 2)*f(iz+4)
...@@ -122,47 +122,47 @@ SUBROUTINE grad_z4(local_Nz,Ngz,inv_deltaz,f,ddz4f) ...@@ -122,47 +122,47 @@ SUBROUTINE grad_z4(local_Nz,Ngz,inv_deltaz,f,ddz4f)
END SUBROUTINE grad_z4 END SUBROUTINE grad_z4
SUBROUTINE interp_z(target,local_Nz,Ngz,f_in,f_out) SUBROUTINE interp_z(target,local_nz,Ngz,f_in,f_out)
! Function meant to interpolate one field defined on a even/odd z into ! Function meant to interpolate one field defined on a even/odd z into
! the other odd/even z grid. ! the other odd/even z grid.
! If Staggered Grid flag (SG) is false, returns identity ! If Staggered Grid flag (SG) is false, returns identity
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: f_in COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: f_in
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: f_out COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: f_out
SELECT CASE(TARGET) SELECT CASE(TARGET)
CASE(1) ! output on even grid CASE(1) ! output on even grid
CALL interp_o2e_z(local_Nz,Ngz,f_in,f_out) CALL interp_o2e_z(local_nz,Ngz,f_in,f_out)
CASE(2) ! output on odd grid CASE(2) ! output on odd grid
CALL interp_e2o_z(local_Nz,Ngz,f_in,f_out) CALL interp_e2o_z(local_nz,Ngz,f_in,f_out)
CASE DEFAULT ! No staggered grid -> usual centered finite differences CASE DEFAULT ! No staggered grid -> usual centered finite differences
f_out = f_in f_out = f_in
END SELECT END SELECT
CONTAINS CONTAINS
SUBROUTINE interp_o2e_z(local_Nz, Ngz,fo,fe) SUBROUTINE interp_o2e_z(local_nz, Ngz,fo,fe)
! gives the value of a field from the odd grid to the even one ! gives the value of a field from the odd grid to the even one
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fo COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fo
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: fe COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fe
INTEGER :: iz INTEGER :: iz
! 4th order interp ! 4th order interp
DO iz = 1,local_Nz DO iz = 1,local_nz
fe(iz) = iz_o2e(-2)*fo(iz ) + iz_o2e(-1)*fo(iz+1) & fe(iz) = iz_o2e(-2)*fo(iz ) + iz_o2e(-1)*fo(iz+1) &
+ iz_o2e( 0)*fo(iz+2) + iz_o2e( 1)*fo(iz+3) + iz_o2e( 0)*fo(iz+2) + iz_o2e( 1)*fo(iz+3)
ENDDO ENDDO
END SUBROUTINE interp_o2e_z END SUBROUTINE interp_o2e_z
SUBROUTINE interp_e2o_z(local_Nz, Ngz,fe,fo) SUBROUTINE interp_e2o_z(local_nz, Ngz,fe,fo)
! gives the value of a field from the even grid to the odd one ! gives the value of a field from the even grid to the odd one
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz, Ngz INTEGER, INTENT(IN) :: local_nz, Ngz
COMPLEX(dp),dimension(local_Nz+Ngz), INTENT(IN) :: fe COMPLEX(dp),dimension(local_nz+Ngz), INTENT(IN) :: fe
COMPLEX(dp),dimension(local_Nz), INTENT(OUT) :: fo COMPLEX(dp),dimension(local_nz), INTENT(OUT) :: fo
INTEGER :: iz INTEGER :: iz
! 4th order interp ! 4th order interp
DO iz = 1,local_Nz DO iz = 1,local_nz
fo(iz) = iz_e2o(-1)*fe(iz+1) + iz_e2o( 0)*fe(iz+2) & fo(iz) = iz_e2o(-1)*fe(iz+1) + iz_e2o( 0)*fe(iz+2) &
+ iz_e2o( 1)*fe(iz+3) + iz_e2o( 2)*fe(iz+4) + iz_e2o( 1)*fe(iz+3) + iz_e2o( 2)*fe(iz+4)
ENDDO ENDDO
...@@ -170,23 +170,23 @@ CONTAINS ...@@ -170,23 +170,23 @@ CONTAINS
END SUBROUTINE interp_z END SUBROUTINE interp_z
SUBROUTINE simpson_rule_z(local_Nz,dz,f,intf) SUBROUTINE simpson_rule_z(local_nz,dz,f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs)) ! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei !from molix BJ Frei
USE prec_const, ONLY: dp, onethird USE prec_const, ONLY: dp, onethird
USE parallel, ONLY: num_procs_z, rank_z, comm_z, manual_0D_bcast USE parallel, ONLY: num_procs_z, rank_z, comm_z, manual_0D_bcast
USE mpi USE mpi
implicit none implicit none
INTEGER, INTENT(IN) :: local_Nz INTEGER, INTENT(IN) :: local_nz
REAL(dp),INTENT(IN) :: dz REAL(dp),INTENT(IN) :: dz
complex(dp),dimension(local_Nz), intent(in) :: f complex(dp),dimension(local_nz), intent(in) :: f
COMPLEX(dp), intent(out) :: intf COMPLEX(dp), intent(out) :: intf
COMPLEX(dp) :: buffer, local_int COMPLEX(dp) :: buffer, local_int
INTEGER :: root, i_, iz, ierr INTEGER :: root, i_, iz, ierr
! Buil local sum using the weights of composite Simpson's rule ! Buil local sum using the weights of composite Simpson's rule
local_int = 0._dp local_int = 0._dp
DO iz = 1,local_Nz DO iz = 1,local_nz
IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz
local_int = local_int + 2._dp*onethird*dz*f(iz) local_int = local_int + 2._dp*onethird*dz*f(iz)
ELSE ! even iz ELSE ! even iz
......
src/control.F90
+ 3
− 2
View file @ a0bebee6
...@@ -61,6 +61,9 @@ SUBROUTINE control ...@@ -61,6 +61,9 @@ SUBROUTINE control
! 2. Main loop ! 2. Main loop
DO DO
CALL cpu_time(t0_step) ! Measuring time CALL cpu_time(t0_step) ! Measuring time
CALL tesend
IF( nlend ) EXIT ! exit do loop
CALL increase_step CALL increase_step
CALL increase_cstep CALL increase_cstep
...@@ -68,8 +71,6 @@ SUBROUTINE control ...@@ -68,8 +71,6 @@ SUBROUTINE control
CALL increase_time CALL increase_time
CALL tesend
IF( nlend ) EXIT ! exit do loop
CALL diagnose(step) CALL diagnose(step)
......
src/diagnose.F90
+ 13
− 18
View file @ a0bebee6
SUBROUTINE diagnose(kstep) SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk ! Diagnostics, writing simulation state to disk
USE basic, ONLY: t0_diag,t1_diag,tc_diag, lu_in, finish, start, cstep, dt, time, tmax, display_h_min_s USE basic, ONLY: t0_diag,t1_diag,tc_diag, lu_in, finish, start, cstep, dt, time, tmax, display_h_min_s
USE diagnostics_par, ONLY: input_fname USE diagnostics_par, ONLY: input_fname
USE processing, ONLY: pflux_x, hflux_x USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id USE parallel, ONLY: my_id
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in) :: kstep INTEGER, INTENT(in) :: kstep
CALL cpu_time(t0_diag) ! Measuring time CALL cpu_time(t0_diag) ! Measuring time
...@@ -24,7 +24,7 @@ SUBROUTINE diagnose(kstep) ...@@ -24,7 +24,7 @@ SUBROUTINE diagnose(kstep)
!! Specific diagnostic calls !! Specific diagnostic calls
CALL diagnose_full(kstep) CALL diagnose_full(kstep)
! Terminal info ! Terminal info
IF ((kstep .GT. 0) .AND. (MOD(cstep, INT(1.0/dt)) == 0) .AND. (my_id .EQ. 0)) THEN IF ((kstep .GE. 0) .AND. (MOD(cstep, INT(1.0/dt)) == 0) .AND. (my_id .EQ. 0)) THEN
WRITE(*,"(A,F6.0,A1,F6.0,A8,G10.2,A8,G10.2,A)")'|t/tmax = ', time,"/",tmax,'| Gxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|' WRITE(*,"(A,F6.0,A1,F6.0,A8,G10.2,A8,G10.2,A)")'|t/tmax = ', time,"/",tmax,'| Gxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|'
ENDIF ENDIF
CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
...@@ -155,7 +155,7 @@ SUBROUTINE diagnose_full(kstep) ...@@ -155,7 +155,7 @@ SUBROUTINE diagnose_full(kstep)
CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 3/)) CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 3/))
CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 2, 4/)) CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,1), (/1, 1, 2, 4/))
! var0d group (gyro transport) ! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN IF (nsave_0d .GT. 0) THEN
CALL creatg(fidres, "/data/var0d", "0d profiles") CALL creatg(fidres, "/data/var0d", "0d profiles")
...@@ -225,13 +225,10 @@ SUBROUTINE diagnose_full(kstep) ...@@ -225,13 +225,10 @@ SUBROUTINE diagnose_full(kstep)
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
END IF END IF
ENDIF ENDIF
!_____________________________________________________________________________ !_____________________________________________________________________________
! 2. Periodic diagnostics ! 2. Periodic diagnostics
! !
IF (kstep .GE. 0) THEN IF (kstep .GE. 0) THEN
! 2.1 0d history arrays ! 2.1 0d history arrays
IF (nsave_0d .GT. 0) THEN IF (nsave_0d .GT. 0) THEN
IF ( MOD(cstep, nsave_0d) == 0 ) THEN IF ( MOD(cstep, nsave_0d) == 0 ) THEN
...@@ -248,15 +245,13 @@ SUBROUTINE diagnose_full(kstep) ...@@ -248,15 +245,13 @@ SUBROUTINE diagnose_full(kstep)
ENDIF ENDIF
ENDIF ENDIF
! 2.4 5d profiles ! 2.4 5d profiles
IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN IF (nsave_5d .GT. 0) THEN
IF (MOD(cstep, nsave_5d) == 0) THEN IF (MOD(cstep, nsave_5d) == 0) THEN
CALL diagnose_5d CALL diagnose_5d
END IF END IF
END IF END IF
!_____________________________________________________________________________ !_____________________________________________________________________________
! 3. Final diagnostics ! 3. Final diagnostics
ELSEIF (kstep .EQ. -1) THEN ELSEIF (kstep .EQ. -1) THEN
CALL attach(fidres, "/data/input","cpu_time",finish-start) CALL attach(fidres, "/data/input","cpu_time",finish-start)
...@@ -446,8 +441,8 @@ SUBROUTINE diagnose_5d ...@@ -446,8 +441,8 @@ SUBROUTINE diagnose_5d
USE fields, ONLY: moments USE fields, ONLY: moments
USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, & USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
local_np, local_nj, local_nky, local_nkx, local_nz, & local_np, local_nj, local_nky, local_nkx, local_nz, &
ngp, ngj, ngz ngp, ngj, ngz, total_na
USE time_integration, ONLY: updatetlevel USE time_integration, ONLY: updatetlevel, ntimelevel
USE diagnostics_par USE diagnostics_par
USE prec_const, ONLY: dp USE prec_const, ONLY: dp
IMPLICIT NONE IMPLICIT NONE
...@@ -470,19 +465,19 @@ SUBROUTINE diagnose_5d ...@@ -470,19 +465,19 @@ SUBROUTINE diagnose_5d
USE parallel, ONLY: gather_pjxyz, num_procs USE parallel, ONLY: gather_pjxyz, num_procs
USE prec_const, ONLY: dp USE prec_const, ONLY: dp
IMPLICIT NONE IMPLICIT NONE
COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: field COMPLEX(dp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text CHARACTER(*), INTENT(IN) :: text
COMPLEX(dp), DIMENSION(local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub COMPLEX(dp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
COMPLEX(dp), DIMENSION(total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full COMPLEX(dp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name CHARACTER(LEN=50) :: dset_name
field_sub = field(1,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),& field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
1:local_nky,1:local_nkx, (1+ngz/2):(local_nz+ngz/2),updatetlevel) 1:local_nky,1:local_nkx, (1+ngz/2):(local_nz+ngz/2),updatetlevel)
field_full = 0; field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN IF (num_procs .EQ. 1) THEN
CALL putarr(fidres, dset_name, field_sub, ionode=0) CALL putarr(fidres, dset_name, field_sub, ionode=0)
ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv) ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
CALL gather_pjxyz(field_sub,field_full,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz) CALL gather_pjxyz(field_sub,field_full,total_na,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
CALL putarr(fidres, dset_name, field_full, ionode=0) CALL putarr(fidres, dset_name, field_full, ionode=0)
ELSE ELSE
CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/)) CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/))
......
src/geometry_mod.F90
+ 9
− 9
View file @ a0bebee6
...@@ -100,7 +100,7 @@ CONTAINS ...@@ -100,7 +100,7 @@ CONTAINS
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient ! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none implicit none
REAL(dp) :: kx,ky REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(local_Nz+Ngz) :: integrant COMPLEX(dp), DIMENSION(local_nz) :: integrant
real(dp) :: G1,G2,G3,Cx,Cy real(dp) :: G1,G2,G3,Cx,Cy
INTEGER :: eo,iz,iky,ikx INTEGER :: eo,iz,iky,ikx
...@@ -138,7 +138,7 @@ CONTAINS ...@@ -138,7 +138,7 @@ CONTAINS
CALL set_kparray(gxx,gxy,gyy,hatB) CALL set_kparray(gxx,gxy,gyy,hatB)
DO eo = 1,Nzgrid DO eo = 1,Nzgrid
! Curvature operator (Frei et al. 2022 eq 2.15) ! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = 1,local_Nz+Ngz DO iz = 1,local_nz+Ngz
G1 = gxx(iz,eo)*gyy(iz,eo)-gxy(iz,eo)*gxy(iz,eo) G1 = gxx(iz,eo)*gyy(iz,eo)-gxy(iz,eo)*gxy(iz,eo)
G2 = gxx(iz,eo)*gyz(iz,eo)-gxy(iz,eo)*gxz(iz,eo) G2 = gxx(iz,eo)*gyz(iz,eo)-gxy(iz,eo)*gxz(iz,eo)
G3 = gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo) G3 = gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo)
...@@ -166,7 +166,7 @@ CONTAINS ...@@ -166,7 +166,7 @@ CONTAINS
CALL set_ikx_zBC_map CALL set_ikx_zBC_map
! !
! Compute the inverse z integrated Jacobian (useful for flux averaging) ! Compute the inverse z integrated Jacobian (useful for flux averaging)
integrant = Jacobian(:,1) ! Convert into complex array integrant = Jacobian(1+ngz/2:local_nz+ngz/2,1) ! Convert into complex array
CALL simpson_rule_z(local_nz,deltaz,integrant,iInt_Jacobian) CALL simpson_rule_z(local_nz,deltaz,integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry END SUBROUTINE eval_magnetic_geometry
...@@ -175,7 +175,7 @@ CONTAINS ...@@ -175,7 +175,7 @@ CONTAINS
! !
SUBROUTINE eval_salpha_geometry SUBROUTINE eval_salpha_geometry
USE grid, ONLY : local_Nz,Ngz,zarray,Nzgrid USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
! evaluate s-alpha geometry model ! evaluate s-alpha geometry model
implicit none implicit none
REAL(dp) :: z REAL(dp) :: z
...@@ -183,7 +183,7 @@ CONTAINS ...@@ -183,7 +183,7 @@ CONTAINS
alpha_MHD = 0._dp alpha_MHD = 0._dp
DO eo = 1,Nzgrid DO eo = 1,Nzgrid
DO iz = 1,local_Nz+Ngz DO iz = 1,local_nz+Ngz
z = zarray(iz,eo) z = zarray(iz,eo)
! metric ! metric
...@@ -228,14 +228,14 @@ CONTAINS ...@@ -228,14 +228,14 @@ CONTAINS
! !
SUBROUTINE eval_zpinch_geometry SUBROUTINE eval_zpinch_geometry
USE grid, ONLY : local_Nz,Ngz,zarray,Nzgrid USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
implicit none implicit none
REAL(dp) :: z REAL(dp) :: z
INTEGER :: iz, eo INTEGER :: iz, eo
alpha_MHD = 0._dp alpha_MHD = 0._dp
DO eo = 1,Nzgrid DO eo = 1,Nzgrid
DO iz = 1,local_Nz+Ngz DO iz = 1,local_nz+Ngz
z = zarray(iz,eo) z = zarray(iz,eo)
! metric ! metric
...@@ -278,7 +278,7 @@ CONTAINS ...@@ -278,7 +278,7 @@ CONTAINS
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
! NOT TESTED ! NOT TESTED
subroutine eval_1D_geometry subroutine eval_1D_geometry
USE grid, ONLY : local_Nz,Ngz,zarray, Nzgrid USE grid, ONLY : local_nz,Ngz,zarray, Nzgrid
! evaluate 1D perp geometry model ! evaluate 1D perp geometry model
implicit none implicit none
REAL(dp) :: z REAL(dp) :: z
...@@ -313,7 +313,7 @@ CONTAINS ...@@ -313,7 +313,7 @@ CONTAINS
USE grid, ONLY: local_nky,Nkx, contains_zmin,contains_zmax, Nexc USE grid, ONLY: local_nky,Nkx, contains_zmin,contains_zmax, Nexc
USE prec_const, ONLY: imagu, pi USE prec_const, ONLY: imagu, pi
IMPLICIT NONE IMPLICIT NONE
! REAL :: shift ! REAL(dp) :: shift
INTEGER :: ikx,iky INTEGER :: ikx,iky
ALLOCATE(ikx_zBC_L(local_nky,Nkx)) ALLOCATE(ikx_zBC_L(local_nky,Nkx))
ALLOCATE(ikx_zBC_R(local_nky,Nkx)) ALLOCATE(ikx_zBC_R(local_nky,Nkx))
......
src/grid_mod.F90
+ 15
− 15
View file @ a0bebee6
...@@ -76,7 +76,7 @@ MODULE grid ...@@ -76,7 +76,7 @@ MODULE grid
REAL(dp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax REAL(dp), PUBLIC, PROTECTED :: local_kxmin, local_kxmax
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: local_zmin, local_zmax
! local z weights for computing simpson rule ! local z weights for computing simpson rule
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zweights_SR
! Numerical diffusion scaling ! Numerical diffusion scaling
REAL(dp), PUBLIC, PROTECTED :: diff_p_coeff, diff_j_coeff REAL(dp), PUBLIC, PROTECTED :: diff_p_coeff, diff_j_coeff
REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
...@@ -381,7 +381,7 @@ CONTAINS ...@@ -381,7 +381,7 @@ CONTAINS
CHARACTER(len=*), INTENT(IN) ::LINEARITY CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx, ikxo INTEGER :: ikx, ikxo
REAL :: Lx_adapted REAL(dp):: Lx_adapted
IF(shear .GT. 0) THEN IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..' IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
! mininal size of box in x to respect dkx = 2pi shear dky ! mininal size of box in x to respect dkx = 2pi shear dky
...@@ -424,7 +424,7 @@ CONTAINS ...@@ -424,7 +424,7 @@ CONTAINS
local_kxmax = 0._dp local_kxmax = 0._dp
DO ikx = ikxs,ikxe DO ikx = ikxs,ikxe
ikxo = ikx - local_nkx_offset ikxo = ikx - local_nkx_offset
kxarray(ikxo) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx))) kxarray(ikxo) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
if (ikx .EQ. Nx/2+1) kxarray(ikxo) = -kxarray(ikxo) if (ikx .EQ. Nx/2+1) kxarray(ikxo) = -kxarray(ikxo)
! Finding kx=0 ! Finding kx=0
IF (kxarray(ikxo) .EQ. 0) THEN IF (kxarray(ikxo) .EQ. 0) THEN
...@@ -441,7 +441,7 @@ CONTAINS ...@@ -441,7 +441,7 @@ CONTAINS
! Build the full grids on process 0 to diagnose it without comm ! Build the full grids on process 0 to diagnose it without comm
! kx ! kx
DO ikx = 1,Nkx DO ikx = 1,Nkx
kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx))) kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx) IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO END DO
ELSE ! Odd number of kx (-2 -1 0 1 2) ELSE ! Odd number of kx (-2 -1 0 1 2)
...@@ -503,8 +503,8 @@ CONTAINS ...@@ -503,8 +503,8 @@ CONTAINS
USE prec_const USE prec_const
USE parallel, ONLY: num_procs_z, rank_z USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE IMPLICIT NONE
REAL :: grid_shift, Lz, zmax, zmin REAL(dp):: grid_shift, Lz, zmax, zmin
INTEGER :: istart, iend, in, Npol, iz, ig, eo, izo INTEGER :: istart, iend, in, Npol, iz, ig, eo
total_nz = Nz total_nz = Nz
! Length of the flux tube (in ballooning angle) ! Length of the flux tube (in ballooning angle)
Lz = 2_dp*pi*Npol Lz = 2_dp*pi*Npol
...@@ -562,10 +562,9 @@ CONTAINS ...@@ -562,10 +562,9 @@ CONTAINS
! Local z array ! Local z array
ALLOCATE(zarray(local_nz+Ngz,Nzgrid)) ALLOCATE(zarray(local_nz+Ngz,Nzgrid))
!! interior point loop !! interior point loop
DO iz = izs,ize DO iz = 1,total_nz
izo = iz - local_nz_offset
DO eo = 1,Nzgrid DO eo = 1,Nzgrid
zarray(izo,eo) = zarray_full(iz) + (eo-1)*grid_shift zarray(iz+ngz/2,eo) = zarray_full(iz) + REAL(eo-1,dp)*grid_shift
ENDDO ENDDO
ENDDO ENDDO
ALLOCATE(local_zmax(Nzgrid),local_zmin(Nzgrid)) ALLOCATE(local_zmax(Nzgrid),local_zmin(Nzgrid))
...@@ -573,16 +572,17 @@ CONTAINS ...@@ -573,16 +572,17 @@ CONTAINS
! Find local extrema ! Find local extrema
local_zmax(eo) = zarray(local_nz+ngz/2,eo) local_zmax(eo) = zarray(local_nz+ngz/2,eo)
local_zmin(eo) = zarray(1+ngz/2,eo) local_zmin(eo) = zarray(1+ngz/2,eo)
print*, zarray
! Fill the ghosts ! Fill the ghosts
DO ig = 1,ngj/2 DO ig = 1,ngz/2
zarray(ig,eo) = local_zmin(eo)-(ngz/2+(ig-1))*deltaz zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),dp)*deltaz
zarray(local_nz+ig,eo) = local_zmax(eo)+ig*deltaz zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,dp)*deltaz
ENDDO ENDDO
! Set up the flags to know if the process contains the tip and/or the tail ! Set up the flags to know if the process contains the tip and/or the tail
! of the z domain (important for z-boundary condition) ! of the z domain (important for z-boundary condition)
IF(abs(local_zmin(eo) - (zmin+(eo-1)*grid_shift)) .LT. EPSILON(zmin)) & IF(abs(local_zmin(eo) - (zmin+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmin)) &
contains_zmin = .TRUE. contains_zmin = .TRUE.
IF(abs(local_zmax(eo) - (zmax+(eo-1)*grid_shift)) .LT. EPSILON(zmax)) & IF(abs(local_zmax(eo) - (zmax+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmax)) &
contains_zmax = .TRUE. contains_zmax = .TRUE.
ENDDO ENDDO
IF(mod(Nz,2) .NE. 0 ) THEN IF(mod(Nz,2) .NE. 0 ) THEN
...@@ -596,7 +596,7 @@ CONTAINS ...@@ -596,7 +596,7 @@ CONTAINS
REAL(dp) :: kx, ky REAL(dp) :: kx, ky
CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2) CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2)
DO eo = 1,Nzgrid DO eo = 1,Nzgrid
DO iz = 1,local_Nz+Ngz DO iz = 1,local_nz+Ngz
DO iky = 1,local_nky DO iky = 1,local_nky
ky = kyarray(iky) ky = kyarray(iky)
DO ikx = 1,local_nkx DO ikx = 1,local_nkx
......
src/inital.F90
+ 3
− 5
View file @ a0bebee6
...@@ -181,7 +181,7 @@ SUBROUTINE init_gyrodens ...@@ -181,7 +181,7 @@ SUBROUTINE init_gyrodens
DO iky=1,local_nky DO iky=1,local_nky
DO iz=1+ngz/2,local_nz+ngz/2 DO iz=1+ngz/2,local_nz+ngz/2
CALL RANDOM_NUMBER(noise) CALL RANDOM_NUMBER(noise)
IF ( (parray(ip) .EQ. 1) .AND. (jarray(ij) .EQ. 1) ) THEN IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp)) moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE ELSE
moments(ia,ip,ij,iky,ikx,iz,:) = 0._dp moments(ia,ip,ij,iky,ikx,iz,:) = 0._dp
...@@ -196,8 +196,6 @@ SUBROUTINE init_gyrodens ...@@ -196,8 +196,6 @@ SUBROUTINE init_gyrodens
ENDIF ENDIF
END DO END DO
END DO END DO
print*, SUM(moments)
stop
! Putting to zero modes that are not in the 2/3 Orszag rule ! Putting to zero modes that are not in the 2/3 Orszag rule
IF (LINEARITY .NE. 'linear') THEN IF (LINEARITY .NE. 'linear') THEN
DO ikx=1,total_nkx DO ikx=1,total_nkx
...@@ -249,7 +247,7 @@ SUBROUTINE init_phi ...@@ -249,7 +247,7 @@ SUBROUTINE init_phi
DO ikx=2,total_nkx/2 DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz) phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO ENDDO
phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz)) !origin must be real phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
END DO END DO
ENDIF ENDIF
!**** ensure no previous moments initialization !**** ensure no previous moments initialization
...@@ -309,7 +307,7 @@ SUBROUTINE init_phi_ppj ...@@ -309,7 +307,7 @@ SUBROUTINE init_phi_ppj
DO ikx=2,total_nkx/2 DO ikx=2,total_nkx/2
phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz) phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
ENDDO ENDDO
phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz)) !origin must be real phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),dp) !origin must be real
END DO END DO
ENDIF ENDIF
!**** ensure no previous moments initialization !**** ensure no previous moments initialization
......
src/memory.F90
+ 26
− 26
View file @ a0bebee6
...@@ -4,7 +4,7 @@ SUBROUTINE memory ...@@ -4,7 +4,7 @@ SUBROUTINE memory
USE array USE array
USE basic, ONLY: allocate_array USE basic, ONLY: allocate_array
USE fields USE fields
USE grid, ONLY: local_Na, local_Np,Ngp ,local_Nj,Ngj, local_nky, local_nkx,local_Nz,Ngz, jmax, pmax USE grid, ONLY: local_Na, local_Np,Ngp ,local_Nj,Ngj, local_nky, local_nkx,local_nz,Ngz, jmax, pmax
USE collision USE collision
USE time_integration, ONLY: ntimelevel USE time_integration, ONLY: ntimelevel
USE prec_const USE prec_const
...@@ -12,30 +12,30 @@ SUBROUTINE memory ...@@ -12,30 +12,30 @@ SUBROUTINE memory
IMPLICIT NONE IMPLICIT NONE
! Electrostatic potential ! Electrostatic potential
CALL allocate_array( phi, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz) CALL allocate_array( phi, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
CALL allocate_array( psi, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz) CALL allocate_array( psi, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
CALL allocate_array(inv_poisson_op, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array(inv_poisson_op, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( inv_ampere_op, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( inv_ampere_op, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( inv_pol_ion, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( inv_pol_ion, 1,local_nky, 1,local_nkx, 1,local_nz)
! CALL allocate_array(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_Nz) ! CALL allocate_array(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_nz)
!Moments related arrays !Moments related arrays
CALL allocate_array( dens, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( dens, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( upar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( upar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( uper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( uper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Tpar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Tpar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Tper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Tper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( temp, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( temp, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Kernel, 1,local_Na, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz, 1,2) CALL allocate_array( Kernel, 1,local_Na, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2)
CALL allocate_array( moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz, 1,ntimelevel ) CALL allocate_array( moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,ntimelevel )
CALL allocate_array( Napjz, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_Nz) CALL allocate_array( Napjz, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nz)
CALL allocate_array( moments_rhs, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz, 1,ntimelevel ) CALL allocate_array( moments_rhs, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz, 1,ntimelevel )
CALL allocate_array( nadiab_moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz+Ngz) CALL allocate_array( nadiab_moments, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz)
CALL allocate_array( ddz_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( ddz_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( ddzND_Napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( ddzND_Napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( interp_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( interp_napj, 1,local_Na, 1,local_Np+Ngp, 1,local_Nj+Ngj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Capj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Capj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Sapj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Sapj, 1,local_Na, 1,local_Np, 1,local_Nj, 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( xnapj, 1,local_Na, 1,local_Np, 1,local_Nj) CALL allocate_array( xnapj, 1,local_Na, 1,local_Np, 1,local_Nj)
CALL allocate_array( xnapp2j, 1,local_Na, 1,local_Np) CALL allocate_array( xnapp2j, 1,local_Na, 1,local_Np)
CALL allocate_array( xnapp1j, 1,local_Na, 1,local_Np) CALL allocate_array( xnapp1j, 1,local_Na, 1,local_Np)
...@@ -67,9 +67,9 @@ SUBROUTINE memory ...@@ -67,9 +67,9 @@ SUBROUTINE memory
!___________________ 2x5D ARRAYS __________________________ !___________________ 2x5D ARRAYS __________________________
!! Collision matrices !! Collision matrices
IF (GK_CO) THEN !GK collision matrices (one for each kperp) IF (GK_CO) THEN !GK collision matrices (one for each kperp)
CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
CALL allocate_array( Caa, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_Nz) CALL allocate_array( Caa, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,local_nky, 1,local_nkx, 1,local_nz)
ELSE !DK collision matrix (same for every k) ELSE !DK collision matrix (same for every k)
CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1) CALL allocate_array( Cab_F, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1)
CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1) CALL allocate_array( Cab_T, 1,Na, 1,Na, 1,(pmax+1)*(jmax+1), 1,(pmax+1)*(jmax+1), 1,1, 1,1, 1,1)
......
src/miller_mod.F90
+ 7
− 7
View file @ a0bebee6
...@@ -6,7 +6,7 @@ MODULE miller ...@@ -6,7 +6,7 @@ MODULE miller
USE basic USE basic
USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0 USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0
! use coordinates,only: gcoor, get_dzprimedz ! use coordinates,only: gcoor, get_dzprimedz
USE grid, ONLY: local_Nky, local_Nkx, local_Nz, Ngz, kyarray, kxarray, zarray, Nz, local_nz_offset USE grid, ONLY: local_Nky, local_Nkx, local_nz, Ngz, kyarray, kxarray, zarray, Nz, local_nz_offset
! use discretization ! use discretization
USE lagrange_interpolation USE lagrange_interpolation
! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol
...@@ -59,12 +59,12 @@ CONTAINS ...@@ -59,12 +59,12 @@ CONTAINS
real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd
real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad. real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad.
real(dp), INTENT(INOUT) :: C_y, C_xy real(dp), INTENT(INOUT) :: C_y, C_xy
real(dp), dimension(1:local_Nz+Ngz,1:2), INTENT(INOUT) :: & real(dp), dimension(1:local_nz+Ngz,1:2), INTENT(INOUT) :: &
gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,& gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
dBdx_,dBdy_,Bfield_,jacobian_,& dBdx_,dBdy_,Bfield_,jacobian_,&
dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, & dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
gradz_coeff_ gradz_coeff_
real(dp), dimension(1:local_Nky,1:local_Nkx,1:local_Nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_ real(dp), dimension(1:local_Nky,1:local_Nkx,1:local_nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
INTEGER :: iz, ikx, iky, eo INTEGER :: iz, ikx, iky, eo
! No parameter in gyacomo yet ! No parameter in gyacomo yet
real(dp) :: sign_Ip_CW=1 ! current sign (only normal current) real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
...@@ -477,7 +477,7 @@ CONTAINS ...@@ -477,7 +477,7 @@ CONTAINS
call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz) call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz)
! Fill the interior of the geom arrays with the results ! Fill the interior of the geom arrays with the results
do eo=1,2 do eo=1,2
DO iz = 1,local_Nz DO iz = 1,local_nz
gxx_(iz+Ngz/2,eo) = gxx_out(iz-local_nz_offset) gxx_(iz+Ngz/2,eo) = gxx_out(iz-local_nz_offset)
gyy_(iz+Ngz/2,eo) = gyy_out(iz-local_nz_offset) gyy_(iz+Ngz/2,eo) = gyy_out(iz-local_nz_offset)
gxz_(iz+Ngz/2,eo) = gxz_out(iz-local_nz_offset) gxz_(iz+Ngz/2,eo) = gxz_out(iz-local_nz_offset)
...@@ -511,7 +511,7 @@ CONTAINS ...@@ -511,7 +511,7 @@ CONTAINS
CALL update_ghosts_z(dxdZ_(:,eo)) CALL update_ghosts_z(dxdZ_(:,eo))
! Curvature operator (Frei et al. 2022 eq 2.15) ! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = 1,local_Nz+Ngz DO iz = 1,local_nz+Ngz
G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo) G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo) G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo) G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
...@@ -536,10 +536,10 @@ CONTAINS ...@@ -536,10 +536,10 @@ CONTAINS
SUBROUTINE update_ghosts_z(fz_) SUBROUTINE update_ghosts_z(fz_)
IMPLICIT NONE IMPLICIT NONE
! INTEGER, INTENT(IN) :: nztot_ ! INTEGER, INTENT(IN) :: nztot_
REAL(dp), DIMENSION(1:local_Nz+Ngz), INTENT(INOUT) :: fz_ REAL(dp), DIMENSION(1:local_nz+Ngz), INTENT(INOUT) :: fz_
REAL(dp), DIMENSION(-2:2) :: buff REAL(dp), DIMENSION(-2:2) :: buff
INTEGER :: status(MPI_STATUS_SIZE), count, last, first INTEGER :: status(MPI_STATUS_SIZE), count, last, first
last = local_Nz+Ngz/2 last = local_nz+Ngz/2
first= 1 + Ngz/2 first= 1 + Ngz/2
IF(Nz .GT. 1) THEN IF(Nz .GT. 1) THEN
IF (num_procs_z .GT. 1) THEN IF (num_procs_z .GT. 1) THEN
......
src/numerics_mod.F90
+ 21
− 20
View file @ a0bebee6
...@@ -72,7 +72,8 @@ END SUBROUTINE build_dv4Hp_table ...@@ -72,7 +72,8 @@ END SUBROUTINE build_dv4Hp_table
SUBROUTINE evaluate_kernels SUBROUTINE evaluate_kernels
USE basic USE basic
USE array, ONLY : kernel!, HF_phi_correction_operator USE array, ONLY : kernel!, HF_phi_correction_operator
USE grid, ONLY : local_Na, local_Nj,Ngj, local_nkx, local_nky, local_nz, Ngz, jarray, kparray USE grid, ONLY : local_Na, local_Nj,Ngj, local_nkx, local_nky, local_nz, Ngz, jarray, kparray,&
nzgrid
USE species, ONLY : sigma2_tau_o2 USE species, ONLY : sigma2_tau_o2
USE prec_const, ONLY: dp USE prec_const, ONLY: dp
IMPLICIT NONE IMPLICIT NONE
...@@ -80,26 +81,26 @@ SUBROUTINE evaluate_kernels ...@@ -80,26 +81,26 @@ SUBROUTINE evaluate_kernels
REAL(dp) :: j_dp, y_, factj REAL(dp) :: j_dp, y_, factj
DO ia = 1,local_Na DO ia = 1,local_Na
DO eo = 1,2 DO eo = 1,nzgrid
DO ikx = 1,local_nkx DO ikx = 1,local_nkx
DO iky = 1,local_nky DO iky = 1,local_nky
DO iz = 1,local_nz + Ngz DO iz = 1,local_nz + Ngz
DO ij = 1, local_Nj + Ngj DO ij = 1, local_nj + Ngj
j_int = jarray(ij) j_int = jarray(ij)
j_dp = REAL(j_int,dp) j_dp = REAL(j_int,dp)
y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2 y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
IF(j_int .LT. 0) THEN IF(j_int .LT. 0) THEN
kernel(ia,ij,iky,ikx,iz,eo) = 0._dp kernel(ia,ij,iky,ikx,iz,eo) = 0._dp
ELSE ELSE
factj = GAMMA(j_dp+1._dp) factj = GAMMA(j_dp+1._dp)
kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
ENDIF ENDIF
ENDDO
! IF (ijs .EQ. 1) & ! if ijs is 1, the ghost kernel has negative index
! kernel(ia,ijgs,iky,ikx,iz,eo) = 0._dp
ENDDO
ENDDO
ENDDO ENDDO
! IF (ijs .EQ. 1) & ! if ijs is 1, the ghost kernel has negative index
! kernel(ia,ijgs,iky,ikx,iz,eo) = 0._dp
ENDDO
ENDDO
ENDDO
ENDDO ENDDO
! !! Correction term for the evaluation of the heat flux ! !! Correction term for the evaluation of the heat flux
! HF_phi_correction_operator(:,:,:) = & ! HF_phi_correction_operator(:,:,:) = &
......
src/parallel_mod.F90
+ 17
− 17
View file @ a0bebee6
...@@ -276,39 +276,39 @@ CONTAINS ...@@ -276,39 +276,39 @@ CONTAINS
!!!!! Gather a field in kinetic + spatial coordinates on rank 0 !!!!! !!!!! Gather a field in kinetic + spatial coordinates on rank 0 !!!!!
!!!!! Gather a field in spatial coordinates on rank 0 !!!!! !!!!! Gather a field in spatial coordinates on rank 0 !!!!!
SUBROUTINE gather_pjxyz(field_loc,field_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot) SUBROUTINE gather_pjxyz(field_loc,field_tot,na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(IN) :: np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot INTEGER, INTENT(IN) :: na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
COMPLEX(dp), DIMENSION(np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc COMPLEX(dp), DIMENSION(na_tot,np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN) :: field_loc
COMPLEX(dp), DIMENSION(np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot COMPLEX(dp), DIMENSION(na_tot,np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p, full y, local z
COMPLEX(dp), DIMENSION(np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p, full y, full z
COMPLEX(dp), DIMENSION(np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_loc) :: buffer_pt_yl_zc ! full p, local y, constant z
COMPLEX(dp), DIMENSION(np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z COMPLEX(dp), DIMENSION(na_tot,np_tot,nky_tot) :: buffer_pt_yt_zc ! full p, full y, constant z
COMPLEX(dp), DIMENSION(np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z COMPLEX(dp), DIMENSION(na_tot,np_loc) :: buffer_pl_cy_zc !local p, constant y, constant z
COMPLEX(dp), DIMENSION(np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z COMPLEX(dp), DIMENSION(na_tot,np_tot) :: buffer_pt_cy_zc ! full p, constant y, constant z
INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij
snd_p = np_loc ! Number of points to send along p (per z,y) snd_p = na_tot*np_loc ! Number of points to send along p (per z,y)
snd_y = np_tot*nky_loc ! Number of points to send along y (per z, full p) snd_y = na_tot*np_tot*nky_loc ! Number of points to send along y (per z, full p)
snd_z = np_tot*nky_tot*nz_loc ! Number of points to send along z (full y, full p) snd_z = na_tot*np_tot*nky_tot*nz_loc ! Number of points to send along z (full y, full p)
root_p = 0; root_z = 0; root_ky = 0 root_p = 0; root_z = 0; root_ky = 0
j: DO ij = 1,nj_tot j: DO ij = 1,nj_tot
x: DO ix = 1,nkx_tot x: DO ix = 1,nkx_tot
z: DO iz = 1,nz_loc z: DO iz = 1,nz_loc
y: DO iy = 1,nky_loc y: DO iy = 1,nky_loc
! fill a buffer to contain a slice of p data at constant j, ky, kx and z ! fill a buffer to contain a slice of p data at constant j, ky, kx and z
buffer_pl_cy_zc(1:np_loc) = field_loc(1:np_loc,ij,iy,ix,iz) buffer_pl_cy_zc(1:na_tot,1:np_loc) = field_loc(1:na_tot,1:np_loc,ij,iy,ix,iz)
CALL MPI_GATHERV(buffer_pl_cy_zc, snd_p, MPI_DOUBLE_COMPLEX, & CALL MPI_GATHERV(buffer_pl_cy_zc, snd_p, MPI_DOUBLE_COMPLEX, &
buffer_pt_cy_zc, rcv_p, dsp_p, MPI_DOUBLE_COMPLEX, & buffer_pt_cy_zc, rcv_p, dsp_p, MPI_DOUBLE_COMPLEX, &
root_p, comm_p, ierr) root_p, comm_p, ierr)
buffer_pt_yl_zc(1:np_tot,iy) = buffer_pt_cy_zc(1:np_tot) buffer_pt_yl_zc(1:na_tot,1:np_tot,iy) = buffer_pt_cy_zc(1:na_tot,1:np_tot)
ENDDO y ENDDO y
! send the full line on p contained by root_p ! send the full line on p contained by root_p
IF(rank_p .EQ. 0) THEN IF(rank_p .EQ. 0) THEN
CALL MPI_GATHERV(buffer_pt_yl_zc, snd_y, MPI_DOUBLE_COMPLEX, & CALL MPI_GATHERV(buffer_pt_yl_zc, snd_y, MPI_DOUBLE_COMPLEX, &
buffer_pt_yt_zc, rcv_yp, dsp_yp, MPI_DOUBLE_COMPLEX, & buffer_pt_yt_zc, rcv_yp, dsp_yp, MPI_DOUBLE_COMPLEX, &
root_ky, comm_ky, ierr) root_ky, comm_ky, ierr)
buffer_pt_yt_zl(1:np_tot,1:nky_tot,iz) = buffer_pt_yt_zc(1:np_tot,1:nky_tot) buffer_pt_yt_zl(1:na_tot,1:np_tot,1:nky_tot,iz) = buffer_pt_yt_zc(1:na_tot,1:np_tot,1:nky_tot)
ENDIF ENDIF
ENDDO z ENDDO z
! send the full line on y contained by root_kyas ! send the full line on y contained by root_kyas
...@@ -319,7 +319,7 @@ CONTAINS ...@@ -319,7 +319,7 @@ CONTAINS
ENDIF ENDIF
! ID 0 (the one who ouptut) rebuild the whole array ! ID 0 (the one who ouptut) rebuild the whole array
IF(my_id .EQ. 0) & IF(my_id .EQ. 0) &
field_tot(1:np_tot,ij,1:nky_tot,ix,1:nz_tot) = buffer_pt_yt_zt(1:np_tot,1:nky_tot,1:nz_tot) field_tot(1:na_tot,1:np_tot,ij,1:nky_tot,ix,1:nz_tot) = buffer_pt_yt_zt(1:na_tot,1:np_tot,1:nky_tot,1:nz_tot)
ENDDO x ENDDO x
ENDDO j ENDDO j
END SUBROUTINE gather_pjxyz END SUBROUTINE gather_pjxyz
......
src/processing_mod.F90
+ 465
− 458
View file @ a0bebee6
This diff is collapsed.
src/solve_EM_fields.F90
+ 3
− 3
View file @ a0bebee6
...@@ -39,7 +39,7 @@ CONTAINS ...@@ -39,7 +39,7 @@ CONTAINS
x:DO ikx = 1,local_nky x:DO ikx = 1,local_nky
y:DO iky = 1,local_nkx y:DO iky = 1,local_nkx
!!!!!!!!!!!!!!! Compute particle charge density q_a n_a for each evolved species !!!!!!!!!!!!!!! Compute particle charge density q_a n_a for each evolved species
DO iz = 1,local_Nz DO iz = 1,local_nz
rho(iz) = 0._dp rho(iz) = 0._dp
DO in=1,total_nj DO in=1,total_nj
DO ia = 1,local_na DO ia = 1,local_na
...@@ -59,7 +59,7 @@ CONTAINS ...@@ -59,7 +59,7 @@ CONTAINS
IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average) IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
! Prepare integrant for z-average ! Prepare integrant for z-average
integrant(iz) = Jacobian(iz+ngz/2,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz) integrant(iz) = Jacobian(iz+ngz/2,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz)
call simpson_rule_z(local_Nz,deltaz,integrant,intf_) ! get the flux averaged phi call simpson_rule_z(local_nz,deltaz,integrant,intf_) ! get the flux averaged phi
fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF ENDIF
rho(iz) = rho(iz) + fsa_phi rho(iz) = rho(iz) + fsa_phi
...@@ -68,7 +68,7 @@ CONTAINS ...@@ -68,7 +68,7 @@ CONTAINS
! IF (ADIAB_I) THEN ! IF (ADIAB_I) THEN
! ENDIF ! ENDIF
!!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!! Inverting the poisson equation
DO iz = 1,local_Nz DO iz = 1,local_nz
phi(iky,ikx,iz+ngz/2) = inv_poisson_op(iky,ikx,iz)*rho(iz) phi(iky,ikx,iz+ngz/2) = inv_poisson_op(iky,ikx,iz)*rho(iz)
ENDDO ENDDO
ENDDO y ENDDO y
......
src/stepon.F90
+ 4
− 3
View file @ a0bebee6
...@@ -6,6 +6,7 @@ SUBROUTINE stepon ...@@ -6,6 +6,7 @@ SUBROUTINE stepon
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
use mpi, ONLY: MPI_COMM_WORLD use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: dp
IMPLICIT NONE IMPLICIT NONE
INTEGER :: num_step, ierr INTEGER :: num_step, ierr
...@@ -140,7 +141,7 @@ SUBROUTINE stepon ...@@ -140,7 +141,7 @@ SUBROUTINE stepon
moments(ia,ip,ij,iky0,ikx,iz,:) = CONJG(moments(ia,ip,ij,iky0,total_nkx+2-ikx,iz,:)) moments(ia,ip,ij,iky0,ikx,iz,:) = CONJG(moments(ia,ip,ij,iky0,total_nkx+2-ikx,iz,:))
END DO END DO
! must be real at origin and top right ! must be real at origin and top right
moments(ia,ip,ij, iky0,ikx0,iz,:) = REAL(moments(ia,ip,ij, iky0,ikx0,iz,:)) moments(ia,ip,ij, iky0,ikx0,iz,:) = REAL(moments(ia,ip,ij, iky0,ikx0,iz,:),dp)
ENDDO a ENDDO a
ENDDO p ENDDO p
ENDDO j ENDDO j
...@@ -150,13 +151,13 @@ SUBROUTINE stepon ...@@ -150,13 +151,13 @@ SUBROUTINE stepon
phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:) phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:)
END DO END DO
! must be real at origin ! must be real at origin
phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:)) phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
! Psi ! Psi
DO ikx=2,total_nkx/2 !symmetry at ky = 0 DO ikx=2,total_nkx/2 !symmetry at ky = 0
psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:) psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:)
END DO END DO
! must be real at origin ! must be real at origin
psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:)) psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
ENDIF ENDIF
END SUBROUTINE enforce_symmetry END SUBROUTINE enforce_symmetry
......
src/tesend.F90
+ 1
− 1
View file @ a0bebee6
...@@ -24,7 +24,7 @@ SUBROUTINE tesend ...@@ -24,7 +24,7 @@ SUBROUTINE tesend
!________________________________________________________________________________ !________________________________________________________________________________
! 2. Test on NRUN ! 2. Test on NRUN
nlend = step .GT. nrun nlend = step .GE. nrun
IF ( nlend ) THEN IF ( nlend ) THEN
CALL speak('NRUN steps done') CALL speak('NRUN steps done')
RETURN RETURN
......
testcases/smallest_problem/fort.90
+ 2
− 2
View file @ a0bebee6
...@@ -91,8 +91,8 @@ ...@@ -91,8 +91,8 @@
&INITIAL_CON &INITIAL_CON
INIT_OPT = 'mom00' INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing' ACT_ON_MODES = 'donothing'
init_background = 1.0 init_background = 0.0
init_noiselvl = 0.0 init_noiselvl = 1.0
iseed = 42 iseed = 42
/ /
&TIME_INTEGRATION_PAR &TIME_INTEGRATION_PAR
......
testcases/smallest_problem/fort_00.90
+ 2
− 2
View file @ a0bebee6
...@@ -73,8 +73,8 @@ ...@@ -73,8 +73,8 @@
&INITIAL_CON &INITIAL_CON
INIT_OPT = 'mom00' INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing' ACT_ON_MODES = 'donothing'
init_background = 1.0 init_background = 0.0
init_noiselvl = 0.0 init_noiselvl = 1.0
iseed = 42 iseed = 42
/ /
&TIME_INTEGRATION_PAR &TIME_INTEGRATION_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