From 3f25a0efc43c0c3ff5407d5519ff62f5fb5c9a10 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 24 May 2022 14:49:34 +0200
Subject: [PATCH] new output start and hyp_z corrections
---
src/array_mod.F90 | 4 ++--
src/auxval.F90 | 22 +++++++++++-----------
src/calculus_mod.F90 | 3 ++-
src/ghosts_mod.F90 | 1 -
src/grid_mod.F90 | 4 +++-
src/memory.F90 | 4 ++--
src/moments_eq_rhs_mod.F90 | 8 ++++----
src/processing_mod.F90 | 10 +++++-----
8 files changed, 29 insertions(+), 27 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 96485a88..75527204 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -14,10 +14,10 @@ MODULE array
! Derivatives and interpolated moments
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz_nepj
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: interp_nepj
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz2_Nepj
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz4_Nepj
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz_nipj
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: interp_nipj
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz2_Nipj
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz4_Nipj
! Arrays to store special initial modes (semi linear simulation)
! Zonal ones (ky=0)
diff --git a/src/auxval.F90 b/src/auxval.F90
index db695e67..34164c59 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -78,16 +78,16 @@ subroutine auxval
ENDDO
CALL mpi_barrier(MPI_COMM_WORLD, ierr)
- IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN
- IF(KIN_E) &
- write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi
- write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi
- ENDIF
- DO ip = ips_i,ipe_i
- DO ij = ijs_i,ije_i
- IF((parray_i(ip)+2*jarray_i(ij) .LE. dmaxi) .AND. (rank_ky + rank_z .EQ. 0))&
- print*, '(',parray_i(ip),',',jarray_i(ij),')'
- ENDDO
- ENDDO
+ ! IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN
+ ! IF(KIN_E) &
+ ! write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi
+ ! write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi
+ ! ENDIF
+ ! DO ip = ips_i,ipe_i
+ ! DO ij = ijs_i,ije_i
+ ! IF((parray_i(ip)+2*jarray_i(ij) .LE. dmaxi) .AND. (rank_ky + rank_z .EQ. 0))&
+ ! print*, '(',parray_i(ip),',',jarray_i(ij),')'
+ ! ENDDO
+ ! ENDDO
END SUBROUTINE auxval
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 76e22a86..466df531 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -16,7 +16,7 @@ MODULE calculus
REAL(dp), dimension(-2:2) :: dz2_usu = &
(/-1.0_dp/12.0_dp, 4.0_dp/3.0_dp, -5.0_dp/2.0_dp, 4.0_dp/3.0_dp, -1.0_dp/12.0_dp /)! 2th derivative, 4th order (for parallel hypdiff)
REAL(dp), dimension(-2:2) :: dz4_usu = &
- (/ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp /)* 0.0625 ! 4th derivative, 2nd order (for parallel hypdiff)
+ (/ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp /) ! 4th derivative, 2nd order (for parallel hypdiff)
PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4
CONTAINS
@@ -99,6 +99,7 @@ SUBROUTINE grad_z4(f,ddz4f)
IF(Nz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
DO iz = izs,ize
ddz4f(iz) = dz4_usu(-2)*f(iz-2) + dz4_usu(-1)*f(iz-1) &
+ +dz4_usu( 0)*f(iz )&
+dz4_usu( 1)*f(iz+1) + dz4_usu( 2)*f(iz+2)
ENDDO
ELSE
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 59cf8ab7..9dff68f4 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -154,7 +154,6 @@ SUBROUTINE update_ghosts_z_e
! last+2 gets first+1
moments_e(:,:,:,:,ize+2,updatetlevel) = moments_e(:,:,:,:,izs+1,updatetlevel)
ENDIF
-
END SUBROUTINE update_ghosts_z_e
SUBROUTINE update_ghosts_z_i
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 36585c36..01d7e1ab 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -405,7 +405,9 @@ CONTAINS
! Z stepping (#interval = #points since periodic)
deltaz = Lz/REAL(Nz,dp)
inv_deltaz = 1._dp/deltaz
- diff_dz_coeff = (deltaz/2._dp)**2
+ diff_dz_coeff = -(deltaz/2._dp)**4 ! adaptive fourth derivative
+ ! non adaptive
+ ! diff_dz_coeff = -1._dp
IF (SG) THEN
grid_shift = deltaz/2._dp
ELSE
diff --git a/src/memory.F90 b/src/memory.F90
index a3837b7c..3b12ca1a 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -33,7 +33,7 @@ SUBROUTINE memory
CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( ddz_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
- CALL allocate_array( ddz2_Nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz4_Nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( interp_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( moments_e_ZF, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, izs,ize)
CALL allocate_array( moments_e_EM, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, izs,ize)
@@ -69,7 +69,7 @@ SUBROUTINE memory
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( ddz_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
- CALL allocate_array( ddz2_Nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz4_Nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( interp_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( moments_i_ZF, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, izs,ize)
CALL allocate_array( moments_i_EM, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, izs,ize)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index d2e316fe..db861edd 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -59,7 +59,7 @@ SUBROUTINE moments_eq_rhs_e
eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
kperp2= kparray(iky,ikx,iz,eo)**2
- IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxe)) THEN
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
!! Compute moments mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
@@ -110,7 +110,7 @@ SUBROUTINE moments_eq_rhs_e
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
- + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) &
+ + mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) &
! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term
@@ -179,7 +179,7 @@ SUBROUTINE moments_eq_rhs_i
p_int = parray_i(ip) ! Hermite degree
eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
kperp2= kparray(iky,ikx,iz,eo)**2
- IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxi)) THEN
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN
!! Compute moments mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
@@ -233,7 +233,7 @@ SUBROUTINE moments_eq_rhs_i
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
- + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) &
+ + mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) &
! Collision term
+ TColl_i(ip,ij,iky,ikx,iz)&
! Nonlinear term
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 5a47ed6f..ffa5ac89 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -166,11 +166,11 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, &
- ddz_nepj, ddz2_Nepj, interp_nepj,&
- ddz_nipj, ddz2_Nipj, interp_nipj
+ ddz_nepj, ddz4_Nepj, interp_nepj,&
+ ddz_nipj, ddz4_Nipj, interp_nipj
USE time_integration, ONLY : updatetlevel
USE model, ONLY : qe_taue, qi_taui, KIN_E, CLOS
- USE calculus, ONLY : grad_z, grad_z2, interp_z
+ USE calculus, ONLY : grad_z, grad_z4, interp_z
IMPLICIT NONE
INTEGER :: eo, p_int, j_int
CALL cpu_time(t0_process)
@@ -239,7 +239,7 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! Compute z derivatives
CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z2(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nepj(ip,ij,iky,ikx,izs:ize))
+ CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz4_Nepj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize))
ENDDO
@@ -257,7 +257,7 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! Compute z first derivative
CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z2(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nipj(ip,ij,iky,ikx,izs:ize))
+ CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz4_Nipj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize))
ENDDO
--
GitLab