From e84a81a4bba706f2c48913c553128a547c5eb941 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 5 May 2023 16:09:16 +0200
Subject: [PATCH] monomial closure done, no convincing results
---
Makefile | 3 +-
fort.90 | 2 +-
src/CLA_mod.F90 | 190 ++++++++++++++++++++++++++++----------------
src/closure_mod.F90 | 98 ++++++++++++++++++-----
src/diagnose.F90 | 6 +-
5 files changed, 207 insertions(+), 92 deletions(-)
diff --git a/Makefile b/Makefile
index 1316ceff..06b6ca82 100644
--- a/Makefile
+++ b/Makefile
@@ -73,7 +73,6 @@ test_svd: F90 = mpif90
test_svd: F90FLAGS = -DTEST_SVD -O3
test_svd: EXEC = $(BINDIR)/gyacomo23_test_svd
test_svd: dirs src/srcinfo.h
-test_svd: dirs src/srcinfo.h
test_svd: compile
# subroutines
dirs:
@@ -148,7 +147,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/closure_mod.o : src/closure_mod.F90 \
$(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o \
- $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
$(OBJDIR)/collision_mod.o : src/collision_mod.F90 \
diff --git a/fort.90 b/fort.90
index b3a5e9fa..88f8dbfb 100644
--- a/fort.90
+++ b/fort.90
@@ -8,7 +8,7 @@
&GRID
pmax = 4
jmax = 2
- Nx = 128
+ Nx = 64
Lx = 200
Ny = 48
Ly = 60
diff --git a/src/CLA_mod.F90 b/src/CLA_mod.F90
index 16acc115..0729a947 100644
--- a/src/CLA_mod.F90
+++ b/src/CLA_mod.F90
@@ -5,8 +5,9 @@ module CLA
! closure
USE prec_const
USE parallel
+ USE FMZM ! For factorial
implicit none
- ! LOCAL VARIABLES
+ ! local variables for DLRA and SVD
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: A_buff,Bf,Br ! buffer and full/reduced rebuilt matrices
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: U ! basis
REAL(xp), DIMENSION(:), ALLOCATABLE :: Sf, Sr ! full and reduced singular values
@@ -15,103 +16,158 @@ module CLA
COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: work
REAL(xp), DIMENSION(:), ALLOCATABLE :: rwork
- REAL(xp), DIMENSION(:), ALLOCATABLE :: ln ! Laguerre coefficients
- REAL(xp), DIMENSION(:), ALLOCATABLE :: hn ! Hermite coefficients
+ ! local variables for monomial truncation (Hermite and Laguerre coeff)
+ ! Coeff for the linear combination of the monomial truncation scheme
+ REAL(xp), PUBLIC, PROTECTED, DIMENSION(:), ALLOCATABLE :: c_Pp1, c_Pp2, c_Jp1
+
#ifdef TEST_SVD
! These routines are meant for testing SVD filtering
PUBLIC :: init_CLA, filter_sv_moments_ky_pj, test_SVD
#endif
- PUBLIC :: invert_utr_matrix, build_hermite_matrix, build_laguerre_matrix
+ PUBLIC :: set_monomial_trunc_coeff
CONTAINS
!--------------Routines to compute coeff for monomial truncation----------
- ! Invert an upper triangular matrix
- SUBROUTINE invert_utr_matrix(U,invU)
+ SUBROUTINE set_monomial_trunc_coeff(Pmax,Jmax)
IMPLICIT NONE
- REAL(xp), DIMENSION(:,:), INTENT(IN) :: U ! Matrix to invert
- REAL(xp), DIMENSION(:,:), INTENT(OUT) :: invU! Result
+ INTEGER, INTENT(IN) :: Pmax,Jmax ! truncation degree (p,j>d are not evolved)
+ REAL(xp), DIMENSION(Pmax+3,Pmax+3) :: HH, iHH ! Hermite matrix to invert (+1 dim for p=0 and +2 dim for coeff P+1,P+2)
+ REAL(xp), DIMENSION(Jmax+2,Jmax+2) :: LL, iLL ! Laguerre matrix to invert (+1 dim for j=0 and +1 dim for coeff J+1)
+ INTEGER :: i,n
+ ! Hermite
+ HH = 0._dp
+ DO i = 1,Pmax+3
+ n = i-1 !poly. degree
+ HH(1:i,i) = get_hermite_coeffs(n)
+ ENDDO
+ CALL inverse_triangular(Pmax+3,HH,iHH)
+
+ ! Laguerre
+ LL = 0._dp
+ DO i = 1,Jmax+2
+ n = i-1 !poly. degree
+ LL(1:i,i) = get_laguerre_coeffs(n)
+ ENDDO
+ CALL inverse_triangular(Jmax+2,LL,iLL)
+
+ ! Allocate and fill the monomial truncation coefficients
+ ALLOCATE(c_Pp1(Pmax+1))
+ DO i = 1,Pmax+1
+ c_Pp1(i) = -iHH(i,Pmax+2)/iHH(Pmax+2,Pmax+2)
+ ENDDO
+ ALLOCATE(c_Pp2(Pmax+1))
+ DO i = 1,Pmax+1
+ c_Pp2(i) = -iHH(i,Pmax+3)/iHH(Pmax+3,Pmax+3)
+ ENDDO
+ ALLOCATE(c_Jp1(Jmax+1))
+ DO i = 1,Jmax+1
+ c_Jp1(i) = 0*iLL(i,Jmax+2)/iLL(Jmax+2,Jmax+2)
+ ENDDO
END SUBROUTINE
- SUBROUTINE build_hermite_matrix(d,U)
+ ! Invert an upper triangular matrix
+ SUBROUTINE inverse_triangular(N,U,invU)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: d
- REAL(xp), DIMENSION(d,d), INTENT(OUT) :: U
+ INTEGER, INTENT(in) :: N ! size of the matrix
+ REAL(xp), DIMENSION(N,N), INTENT(IN) :: U ! Matrix to invert
+ REAL(xp), DIMENSION(N,N), INTENT(OUT) :: invU! Result
+ ! local variables
+ INTEGER :: info
+ invU = U
+#ifdef SINGLE_PRECISION
+ CALL STRTRI('U','N',N,invU,N,info)
+#else
+ CALL DTRTRI('U','N',N,invU,N,info)
+#endif
+ IF (info .LT. 0) THEN
+ print*, info
+ ERROR STOP "STOP, LAPACK TRTRI: illegal value at index"
+ ELSEIF(info .GT. 0) THEN
+ ERROR STOP "STOP, LAPACK TRTRI: singular matrix"
+ ENDIF
END SUBROUTINE
- SUBROUTINE build_laguerre_matrix(d,U)
+ ! Compute the kth coefficient of the nth degree Hermite
+ ! adapted from LaguerrePoly.m by David Terr, Raytheon, 5-10-04
+ FUNCTION get_hermite_coeffs(n) RESULT(hn)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: d
- REAL(xp), DIMENSION(d,d), INTENT(OUT) :: U
-
-
- END SUBROUTINE
+ INTEGER, INTENT(IN) :: n ! Polyn. Deg. and index
+ REAL(xp), DIMENSION(n+1) :: hn ! Hermite coefficients
+ REAL(xp), DIMENSION(n+1) :: hnm2, hnm1, tmp
+ INTEGER :: k, e, factn
+ IF (n .EQ. 0) THEN
+ hn = 1._xp
+ ELSEIF (n .EQ. 1) THEN
+ hn = [2._xp, 0._xp]
+ ELSE
+ hnm2 = 0._xp
+ hnm2(n+1) = 1._xp
+ hnm1 = 0._xp
+ hnm1(n) = 2._xp
+ DO k = 2,n
+ hn = 0._xp
+ DO e = n-k+1,n,2
+ hn(e) = REAL(2*(hnm1(e+1) - (k-1)*hnm2(e)),xp)
+ END DO
+ hn(n+1) = REAL(-2*(k-1)*hnm2(n+1),xp)
+ IF (k .LT. n) THEN
+ hnm2 = hnm1
+ hnm1 = hn
+ END IF
+ END DO
+ END IF
+ ! Normalize the polynomials
+ factn = 1
+ DO k = 1,n
+ factn = k*factn
+ ENDDO
+ hn = hn/SQRT(REAL(2**n*factn,xp))
+ ! reverse order
+ tmp = hn
+ DO k = 1,n+1
+ hn(n+1-(k-1)) = tmp(k)
+ ENDDO
+ END FUNCTION
! Compute the kth coefficient of the nth degree Laguerre
! adapted from LaguerrePoly.m by David Terr, Raytheon, 5-11-04
! lnk = (-1)^k/k! binom(n,k), s.t. Ln(x) = sum_k lnk x^k
- SUBROUTINE set_laguerre_coeffs(n)
+ FUNCTION get_laguerre_coeffs(n) RESULT(ln)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: n ! Polyn. Deg. and index
- INTEGER, DIMENSION(n+1) :: lnm2, lnm1
+ INTEGER, INTENT(IN) :: n ! Polyn. Deg. and index
+ REAL(xp), DIMENSION(n+1) :: ln ! Laguerre coefficients
+ REAL(xp), DIMENSION(n+1) :: lnm2, lnm1, tmp
INTEGER :: k, e
- ALLOCATE(ln(n+1))
- IF (n == 0) THEN
- ln = 1
- ELSEIF (n == 1) THEN
- ln = [-1, 1]
+ IF (n .EQ. 0) THEN
+ ln = 1._xp
+ ELSEIF (n .EQ. 1) THEN
+ ln = [-1._xp, 1._xp]
ELSE
- lnm2 = 0.0
- lnm2(n+1) = 1
- lnm1 = 0.0
- lnm1(n) = -1
- lnm1(n+1) = 1
+ lnm2 = 0._xp
+ lnm2(n+1) = 1._xp
+ lnm1 = 0._xp
+ lnm1(n) =-1._xp
+ lnm1(n+1) = 1._xp
DO k = 2, n
- ln = 0.0
+ ln = 0
DO e = n-k+1, n
- ln(e) = (2*k-1)*lnm1(e) - lnm1(e+1) + (1-k)*lnm2(e)
+ ln(e) = REAL((2*k-1)*lnm1(e) - lnm1(e+1) + (1-k)*lnm2(e),xp)
END DO
- ln(n+1) = (2*k-1)*lnm1(n+1) + (1-k)*lnm2(n+1)
- ln = ln/k
+ ln(n+1) = REAL((2*k-1)*lnm1(n+1) + (1-k)*lnm2(n+1),xp)
+ ln = ln/REAL(k,xp)
IF (k < n) THEN
lnm2 = lnm1
lnm1 = ln
END IF
END DO
END IF
- END SUBROUTINE
-
- ! Compute the kth coefficient of the nth degree hermite
- ! (-1)^k/k! binom(j,k)
- SUBROUTINE set_hermite_coeffs(n)
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: n ! Polyn. Deg. and index
- INTEGER, DIMENSION(n+1) :: hnm2, hnm1
- INTEGER :: k, e
- ALLOCATE(hn(n+1))
- IF (n .EQ. 0) THEN
- hn = 1
- ELSEIF (n .EQ. 1) THEN
- hn = [2, 0]
- ELSE
- hnm2 = 0
- hnm2(n+1) = 1
- hnm1 = 0
- hnm1(n) = 2
-
- DO k = 2, n
- hn = 0
- DO e = n-k+1, 2, -2
- hn(e) = 2*(hnm1(e+1) - (k-1)*hnm2(e))
- END DO
- hn(n+1) = -2*(k-1)*hnm2(n+1)
- IF (k < n) THEN
- hnm2 = hnm1
- hnm1 = hn
- END IF
- END DO
- END IF
- END SUBROUTINE
+ ! reverse order
+ tmp = ln
+ DO k = 1,n+1
+ ln(n+1-(k-1)) = tmp(k)
+ ENDDO
+ END FUNCTION
!--------------Routines to test singular value filtering -----------------
#ifdef TEST_SVD
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index a27cd910..299fc626 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -23,8 +23,9 @@ END SUBROUTINE
SUBROUTINE set_closure_model
USE basic, ONLY: speak
- USE grid, ONLY: local_np, ngp, local_nj, ngj, parray, jarray,&
- pmax, jmax
+ USE grid, ONLY: local_np, ngp, local_nj, ngj, parray, jarray,&
+ pmax, jmax
+ USE CLA, ONLY: set_monomial_trunc_coeff
IMPLICIT NONE
INTEGER :: ip,ij
! adapt the dmax if it is set <0
@@ -36,7 +37,14 @@ SUBROUTINE set_closure_model
! set the evolve mom array
ALLOCATE(evolve_mom(local_np+ngp,local_nj+ngj))
SELECT CASE(hierarchy_closure)
- CASE('truncation','monomial')
+ CASE('truncation')
+ DO ip = 1,local_np+ngp
+ DO ij = 1, local_nj+ngj
+ evolve_mom(ip,ij) = ((parray(ip).GE.0) .AND. (jarray(ij).GE.0)) &
+ .AND. (parray(ip).LE.pmax) .AND. (jarray(ij).LE.jmax)
+ ENDDO
+ ENDDO
+ CASE('monomial')
DO ip = 1,local_np+ngp
DO ij = 1, local_nj+ngj
evolve_mom(ip,ij) = ((parray(ip).GE.0) .AND. (jarray(ij).GE.0)) &
@@ -51,8 +59,14 @@ SUBROUTINE set_closure_model
ENDDO
ENDDO
CASE DEFAULT
- ERROR STOP "closure scheme not recognized (avail: truncation,max_degree)"
+ ERROR STOP "closure scheme not recognized (avail: truncation,max_degree,monomial)"
END SELECT
+ ! If monomial truncation, setup the coefficients required
+ SELECT CASE(hierarchy_closure)
+ CASE('monomial')
+ CALL set_monomial_trunc_coeff(pmax,jmax)
+ END SELECT
+
! Set the nonlinear closure scheme (truncation of sum over n in Sapj)
ALLOCATE(nmaxarray(local_nj))
SELECT CASE(nonlinear_closure)
@@ -74,31 +88,79 @@ SUBROUTINE set_closure_model
CASE DEFAULT
ERROR STOP "nonlinear closure scheme not recognized (avail: truncation,anti_laguerre_aliasing,full_sum)"
END SELECT
-
- ! If monomial truncation, setup the coefficients required
- SELECT CASE(nonlinear_closure)
- CASE('monomial')
- END SELECT
-
END SUBROUTINE set_closure_model
! Positive Oob indices are approximated with a model
SUBROUTINE apply_closure_model
USE prec_const, ONLY: xp
- USE grid, ONLY: local_nj,ngj,local_np,ngp,local_na
+ USE grid, ONLY: local_nj,ngj,local_np,ngp,local_na,deltap
USE fields, ONLY: moments
USE time_integration, ONLY: updatetlevel
+ USE CLA, ONLY: c_Pp1, c_Pp2, c_Jp1
+ USE parallel, ONLY: num_procs_p
IMPLICIT NONE
- INTEGER ::ij,ip,ia
+ INTEGER ::ij,ip,ia, ipi, iji, ipmax, ijmax, iq, ik, iqi, iki
+ ! Set to zero all non evolving moments
+ DO ij = 1, local_nj+ngj
+ DO ip = 1,local_np+ngp
+ DO ia = 1,local_na
+ IF(.NOT. evolve_mom(ip,ij))&
+ moments(ia,ip,ij,:,:,:,updatetlevel) = 0._xp
+ ENDDO
+ ENDDO
+ ENDDO
SELECT CASE (hierarchy_closure)
CASE('truncation','max_degree')
- DO ij = 1, local_nj+ngj
- DO ip = 1,local_np+ngp
- DO ia = 1,local_na
- IF(.NOT. evolve_mom(ip,ij))&
- moments(ia,ip,ij,:,:,:,updatetlevel) = 0._xp
+ ! do nothing
+ CASE('monomial')
+ IF(num_procs_p .GT. 1) ERROR STOP "STOP: monomial closure is not parallelized in p"
+ !
+ ipmax = local_np+ngp/2
+ ijmax = local_nj+ngj/2
+ DO ia = 1,local_na
+ IF(deltap .EQ. 1) THEN ! We solve for every Hermite
+ !! P+1 truncation : NaP+1,j = sum_{q=0}^P c^{P+1}_q Nqj
+ DO ij = 1, local_nj
+ iji = ij + ngj/2
+ moments(ia,ipmax+1,iji,:,:,:,updatetlevel) = 0._xp
+ DO iq = 1, local_np
+ iqi = iq + ngp/2
+ moments(ia,ipmax+1,iji,:,:,:,updatetlevel) = moments(ia,ipmax+1,iji,:,:,:,updatetlevel) &
+ + c_Pp1(iq) * moments(ia,iqi,iji,:,:,:,updatetlevel)
+ ENDDO
+ ENDDO
+ !! P+2 truncation : NaP+2,j = sum_{q=0}^P c^{P+2}_q Nqj
+ DO ij = 1, local_nj
+ iji = ij + ngj/2
+ moments(ia,ipmax+2,iji,:,:,:,updatetlevel) = 0._xp
+ DO iq = 1, local_np
+ iqi = iq + ngp/2
+ moments(ia,ipmax+2,iji,:,:,:,updatetlevel) = moments(ia,ipmax+2,iji,:,:,:,updatetlevel) &
+ + c_Pp2(iq) * moments(ia,iqi,iji,:,:,:,updatetlevel)
+ ENDDO
+ ENDDO
+ ELSE
+ !! P+2 truncation only : NaP+2,j = sum_{q=0}^P c^{P+2}_q Nqj
+ DO ij = 1, local_nj
+ iji = ij + ngj/2
+ moments(ia,ipmax+1,iji,:,:,:,updatetlevel) = 0._xp
+ DO iq = 1, local_np
+ iqi = iq + ngp/2
+ moments(ia,ipmax+1,iji,:,:,:,updatetlevel) = moments(ia,ipmax+1,iji,:,:,:,updatetlevel) &
+ + c_Pp2(2*(iq-1)+1) * moments(ia,iqi,iji,:,:,:,updatetlevel)
+ ENDDO
+ ENDDO
+ ENDIF
+ !! J+1 truncation : Nap,J+1 = sum_{k=0}^J c^{J+1}_k Npk
+ DO ip = 1, local_np
+ ipi = ip + ngp/2
+ moments(ia,ipi,ijmax+1,:,:,:,updatetlevel) = 0._xp
+ DO ik = 1, local_nj
+ iki = ik + ngj/2
+ moments(ia,iki,ijmax+1,:,:,:,updatetlevel) = moments(ia,iki,ijmax+1,:,:,:,updatetlevel) &
+ + c_Jp1(ik) * moments(ia,ipi,iki,:,:,:,updatetlevel)
ENDDO
- ENDDO
+ ENDDO
ENDDO
CASE DEFAULT
ERROR STOP "closure scheme not recognized"
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index f9e3aa8d..d5aca960 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -96,7 +96,6 @@ SUBROUTINE diagnose_full(kstep)
USE basic, ONLY: speak,chrono_runt,&
cstep,iframe0d,iframe3d,iframe5d,crashed
USE grid, ONLY: &
- local_nj,local_nky,local_nkx,local_nz,ngj,ngz,&
parray_full,pmax,jarray_full,jmax,&
kyarray_full,kxarray_full,zarray_full
USE diagnostics_par
@@ -105,7 +104,6 @@ SUBROUTINE diagnose_full(kstep)
USE model, ONLY: EM
USE parallel, ONLY: my_id, comm0
USE collision, ONLY: coll_outputinputs
- USE geometry, ONLY: gxx,gxy,gyy,gxz,gyz,gzz,hatR,hatZ,hatB,dBdx,dBdy,dBdz,Jacobian,gradz_coeff,Ckxky
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
INTEGER, parameter :: BUFSIZE = 2
@@ -282,7 +280,7 @@ SUBROUTINE diagnose_txtonly(kstep)
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
INTEGER, parameter :: BUFSIZE = 2
- INTEGER :: rank = 0, ierr
+ INTEGER :: rank = 0
INTEGER :: dims(1) = (/0/)
IF (kstep .GE. 0) THEN
! output the transport in a txt file
@@ -507,7 +505,7 @@ SUBROUTINE diagnose_5d
CONTAINS
SUBROUTINE write_field5d(field, text)
- USE basic, ONLY: GATHERV_OUTPUT, jobnum, dt
+ USE basic, ONLY: jobnum, dt
USE futils, ONLY: attach, putarr!, putarrnd
USE parallel, ONLY: gather_pjxyz, num_procs
USE prec_const, ONLY: xp
--
GitLab