monomial closure done, no convincing results

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 \

fort.90

@@ -8,7 +8,7 @@
 &GRID
   pmax   = 4
   jmax   = 2
-  Nx     = 128
+  Nx     = 64
   Lx     = 200
   Ny     = 48
   Ly     = 60

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

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"

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