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

GF closure implemented

parent d118c2b7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/advance_field.F90
+ 13
− 5
View file @ ffd78d08
......@@ -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)
......
src/closure_mod.F90
+ 43
− 36
View file @ ffd78d08
......@@ -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
src/compute_Sapj.F90
+ 11
− 20
View file @ ffd78d08
......@@ -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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......
src/moments_eq_rhs.F90
+ 15
− 6
View file @ ffd78d08
......@@ -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
......
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