From ffd78d087fdde8999e330bf5e802453aa448395d Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 8 Oct 2021 14:03:32 +0200
Subject: [PATCH] GF closure implemented
---
src/advance_field.F90 | 18 +++++++---
src/closure_mod.F90 | 79 +++++++++++++++++++++++-------------------
src/compute_Sapj.F90 | 31 ++++++-----------
src/moments_eq_rhs.F90 | 21 +++++++----
4 files changed, 82 insertions(+), 67 deletions(-)
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 9b7f7938..4f0aa5b9 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -20,20 +20,28 @@ CONTAINS
USE time_integration
USE grid
use prec_const
+ USE model, ONLY: CLOS
use fields, ONLY: moments_e, moments_i
use array, ONLY: moments_rhs_e, moments_rhs_i
IMPLICIT NONE
+ INTEGER :: p_int, j_int
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- CALL advance_field(moments_i(ip,ij,:,:,:,:), moments_rhs_i(ip,ij,:,:,:,:))
- ENDDO
- ENDDO
+ CALL cpu_time(t0_adv_field)
DO ip=ips_e,ipe_e
+ p_int = parray_e(ip)
DO ij=ijs_e,ije_e
+ IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
CALL advance_field(moments_e(ip,ij,:,:,:,:), moments_rhs_e(ip,ij,:,:,:,:))
ENDDO
ENDDO
+ DO ip=ips_i,ipe_i
+ p_int = parray_i(ip)
+ DO ij=ijs_i,ije_i
+ j_int = jarray_i(ij)
+ IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxi))&
+ CALL advance_field(moments_i(ip,ij,:,:,:,:), moments_rhs_i(ip,ij,:,:,:,:))
+ ENDDO
+ ENDDO
! Execution time end
CALL cpu_time(t1_adv_field)
tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field)
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 1b4dabc5..c8055420 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -8,7 +8,7 @@ USE fields, ONLY: moments_e, moments_i
USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
-PUBLIC :: apply_closure_model
+PUBLIC :: apply_closure_model, ghosts_truncation
CONTAINS
@@ -17,42 +17,51 @@ SUBROUTINE apply_closure_model
IMPLICIT NONE
CALL cpu_time(t0_clos)
-! zero truncation, An+1=0 for n+1>nmax
IF (CLOS .EQ. 0) THEN
+ ! zero truncation, An+1=0 for n+1>nmax only
+ CALL ghosts_truncation
+
+
+ ELSEIF (CLOS .EQ. 1) THEN
+ ! Truncation at highest fully represented kinetic moment
+ ! e.g. Dmax = 3 means
+ ! all Napj s.t. p+2j <= 3
+ ! -> (p,j) allowed are (0,0),(1,0),(0,1),(2,0),(1,1),(3,0)
+ ! =>> Dmax is Pmax, condition is p+2j<=Pmax
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
DO ip = ipsg_e,ipeg_e
- moments_e(ip,ijsg_e,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp
- ENDDO
- DO ij = ijsg_e,ijeg_e
- moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipsg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ DO ij = ijsg_e,ijeg_e
+ IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) &
+ moments_e(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ ENDDO
ENDDO
- kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
- kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
-
DO ip = ipsg_i,ipeg_i
- moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp
+ DO ij = ijsg_i,ijeg_i
+ IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) &
+ moments_i(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ ENDDO
ENDDO
- DO ij = ijsg_i,ijeg_i
- moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipsg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- ENDDO
- kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
- kernel_i(ijeg_i,ikx,iky,iz) = 0._dp
ENDDO
ENDDO
ENDDO
+ ! + ghosts truncation
+ CALL ghosts_truncation
+ ELSE
+ if(my_id .EQ. 0) write(*,*) '! Closure scheme not found !'
- ! zero truncation, An+1=0 for n+1>nmax
- ELSEIF (CLOS .EQ. 1) THEN
+ ENDIF
+
+ CALL cpu_time(t1_clos)
+ tc_clos = tc_clos + (t1_clos - t0_clos)
+END SUBROUTINE apply_closure_model
+
+! Positive Oob indices are approximated with a model
+SUBROUTINE ghosts_truncation
+ IMPLICIT NONE
+
+! zero truncation, An+1=0 for n+1>nmax
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
@@ -61,10 +70,12 @@ SUBROUTINE apply_closure_model
moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
DO ij = ijsg_e,ijeg_e
- moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_e(ipsg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
+ moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ ENDIF
ENDDO
kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
@@ -74,23 +85,19 @@ SUBROUTINE apply_closure_model
moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
DO ij = ijsg_i,ijeg_i
- moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_i(ipsg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ IF(deltapi .EQ. 1) THEN ! Must truncate the second stencil
+ moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ ENDIF
ENDDO
kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
kernel_i(ijeg_i,ikx,iky,iz) = 0._dp
ENDDO
ENDDO
ENDDO
- ELSE
- if(my_id .EQ. 0) write(*,*) '! Closure scheme not found !'
-
- ENDIF
- CALL cpu_time(t1_clos)
- tc_clos = tc_clos + (t1_clos - t0_clos)
-END SUBROUTINE apply_closure_model
+END SUBROUTINE ghosts_truncation
END module closure
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 107e3349..5a342d32 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -19,23 +19,20 @@ SUBROUTINE compute_Sapj
REAL(dp), DIMENSION(ixs:ixe,iys:iye) :: fr_real, gz_real
REAL(dp), DIMENSION(ixs:ixe,iys:iye) :: fz_real, gr_real, f_times_g
- INTEGER :: in, is
+ INTEGER :: in, is, p_int, j_int
INTEGER :: nmax, smax ! Upper bound of the sums
REAL(dp):: kx, ky, kerneln
- LOGICAL :: COMPUTE_ONLY_EVEN_P = .true.
! Execution time start
CALL cpu_time(t0_Sapj)
-! If we have a parallel dynamic, odd p are coupled with even ones
-IF(Nz .GT. 1) COMPUTE_ONLY_EVEN_P = .false.
-
zloop: DO iz = izs,ize
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
-
- ! we check if poly degree is even (eq to index is odd) to spare computation
- IF ((MODULO(ip,2) .EQ. 1) .OR. (.NOT. COMPUTE_ONLY_EVEN_P)) THEN
+ p_int = parray_e(ip)
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
+ j_int=jarray_e(ij)
+ ! GF closure check (spare computations too)
+ GF_CLOSURE_e: IF ((CLOS.EQ.1) .AND. (p_int+2*j_int .GT. dmaxe)) THEN
real_data_c = 0._dp ! initialize sum over real nonlinear term
@@ -112,12 +109,8 @@ zloop: DO iz = izs,ize
Sepj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
+ ENDIF GF_CLOSURE_e
ENDDO jloope
-
- ELSE
- ! Cancel the non lin term if we are dealing with odd Hermite degree
- Sepj(ip,:,:,:,iz) = 0._dp
- ENDIF
ENDDO ploope
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -125,9 +118,11 @@ zloop: DO iz = izs,ize
ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
! we check if poly degree is even (eq to index is odd) to spare computation
- IF (MODULO(ip,2) .EQ. 1 .OR. (.NOT. COMPUTE_ONLY_EVEN_P)) THEN
-
+ !EVEN_P_i: IF (.TRUE. .OR. (MODULO(ip,2) .EQ. 1) .OR. (.NOT. COMPUTE_ONLY_EVEN_P)) THEN
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
+ j_int=jarray_i(ij)
+ ! GF closure check (spare computations too)
+ GF_CLOSURE_i: IF ((CLOS.EQ.1) .AND. (p_int+2*j_int .GT. dmaxi)) THEN
real_data_c = 0._dp ! initialize sum over real nonlinear term
! Set non linear sum truncation
@@ -203,12 +198,8 @@ zloop: DO iz = izs,ize
Sipj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
-
+ ENDIF GF_CLOSURE_i
ENDDO jloopi
- ELSE
- ! Cancel the non lin term if we are dealing with odd Hermite degree
- Sipj(ip,:,:,:,iz) = 0._dp
- ENDIF
ENDDO ploopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 24a28360..ce74ae12 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -31,14 +31,17 @@ SUBROUTINE moments_eq_rhs_e
CALL cpu_time(t0_rhs)
ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices
- p_int= parray_e(ip) ! Hermite polynom. degree
-
+ p_int = parray_e(ip) ! Hermite polynom. degree
delta_p0 = 0._dp; delta_p1 = 0._dp; delta_p2 = 0._dp
IF(p_int .EQ. 0) delta_p0 = 1._dp
IF(p_int .EQ. 1) delta_p1 = 1._dp
IF(p_int .EQ. 2) delta_p2 = 1._dp
jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
+ j_int = jarray_e(ij)
+ GF_CLOSURE_e : IF( (CLOS .EQ. 1) .AND. (p_int+2*j_int .GT. Dmaxe) ) THEN
+ !skip
+ ELSE
! Loop on kspace
zloope : DO iz = izs,ize
! Periodic BC for first order derivative
@@ -60,9 +63,9 @@ SUBROUTINE moments_eq_rhs_e
Tnepj = xnepj(ip,ij) * (moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0)
! term propto n_e^{p+2,j}
- Tnepp2j = xnepp2j(ip) * moments_e(ip+2,ij,ikx,iky,iz,updatetlevel)
+ Tnepp2j = xnepp2j(ip) * moments_e(ip+(2/deltape),ij,ikx,iky,iz,updatetlevel)
! term propto n_e^{p-2,j}
- Tnepm2j = xnepm2j(ip) * (moments_e(ip-2,ij,ikx,iky,iz,updatetlevel) &
+ Tnepm2j = xnepm2j(ip) * (moments_e(ip-(2/deltape),ij,ikx,iky,iz,updatetlevel) &
+kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p2)
! term propto n_e^{p,j+1}
Tnepjp1 = xnepjp1(ij) * (moments_e(ip,ij+1,ikx,iky,iz,updatetlevel) &
@@ -130,6 +133,7 @@ SUBROUTINE moments_eq_rhs_e
END DO kxloope
END DO zloope
+ ENDIF GF_CLOSURE_e
END DO jloope
END DO ploope
@@ -178,6 +182,10 @@ SUBROUTINE moments_eq_rhs_i
IF(p_int .EQ. 2) delta_p2 = 1._dp
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
+ j_int = jarray_i(ij)
+ GF_CLOSURE_i : IF( (CLOS .EQ. 1) .AND. (p_int+2*j_int .GT. Dmaxe) ) THEN
+ !skip
+ ELSE
! Loop on kspace
zloopi : DO iz = izs,ize
! Periodic BC for first order derivative
@@ -199,9 +207,9 @@ SUBROUTINE moments_eq_rhs_i
Tnipj = xnipj(ip,ij) * (moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0)
! term propto n_i^{p+2,j}
- Tnipp2j = xnipp2j(ip) * moments_i(ip+2,ij,ikx,iky,iz,updatetlevel)
+ Tnipp2j = xnipp2j(ip) * moments_i(ip+(2/deltapi),ij,ikx,iky,iz,updatetlevel)
! term propto n_i^{p-2,j}
- Tnipm2j = xnipm2j(ip) * (moments_i(ip-2,ij,ikx,iky,iz,updatetlevel) &
+ Tnipm2j = xnipm2j(ip) * (moments_i(ip-(2/deltapi),ij,ikx,iky,iz,updatetlevel) &
+kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p2)
! term propto n_e^{p,j+1}
Tnipjp1 = xnipjp1(ij) * (moments_i(ip,ij+1,ikx,iky,iz,updatetlevel) &
@@ -269,6 +277,7 @@ SUBROUTINE moments_eq_rhs_i
END DO kxloopi
END DO zloopi
+ ENDIF GF_CLOSURE_i
END DO jloopi
END DO ploopi
--
GitLab