From 48240be2496bcdb55431f3e69b088d23e19780ee Mon Sep 17 00:00:00 2001
From: ahoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 4 May 2023 22:41:07 +0200
Subject: [PATCH] Renaming the DLRA module CLA and start of monomial trunc
---
Makefile | 10 +-
src/{DLRA_mod.F90 => CLA_mod.F90} | 126 ++++++++++++++++--
src/auxval.F90 | 4 +-
src/basic_mod.F90 | 4 +-
src/closure_mod.F90 | 9 +-
src/diagnose.F90 | 6 +-
src/stepon.F90 | 2 +-
testcases/DLRA_zpinch/base_case/fort_00.90 | 2 +-
testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 | 2 +-
testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 | 2 +-
10 files changed, 136 insertions(+), 31 deletions(-)
rename src/{DLRA_mod.F90 => CLA_mod.F90} (71%)
diff --git a/Makefile b/Makefile
index af8f6958..b798962f 100644
--- a/Makefile
+++ b/Makefile
@@ -112,7 +112,7 @@ $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/parallel_mod.o
$(OBJDIR)/ppexit.o $(OBJDIR)/prec_const_mod.o \
$(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/restarts_mod.o \
$(OBJDIR)/solve_EM_fields.o $(OBJDIR)/species_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
-$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
+$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
# To compile the executable
compile: $(FOBJ)
@@ -132,7 +132,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
$(OBJDIR)/auxval.o : src/auxval.F90 \
$(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
$(OBJDIR)/geometry_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o \
- $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/DLRA_mod.o
+ $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
$(OBJDIR)/basic_mod.o : src/basic_mod.F90 \
@@ -302,7 +302,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
$(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o \
$(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o\
$(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
- $(OBJDIR)/DLRA_mod.o
+ $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
$(OBJDIR)/tesend.o : src/tesend.F90 \
@@ -317,7 +317,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
$(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
- $(OBJDIR)/DLRA_mod.o : src/DLRA_mod.F90 \
+ $(OBJDIR)/CLA_mod.o : src/CLA_mod.F90 \
$(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/time_integration_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/DLRA_mod.F90 -o $@
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/CLA_mod.F90 -o $@
diff --git a/src/DLRA_mod.F90 b/src/CLA_mod.F90
similarity index 71%
rename from src/DLRA_mod.F90
rename to src/CLA_mod.F90
index 3be00e96..16acc115 100644
--- a/src/DLRA_mod.F90
+++ b/src/CLA_mod.F90
@@ -1,8 +1,11 @@
-module DLRA
+module CLA
+ ! This si the Computational Linear Algebra module.
+ ! It contains tools as routines to DO singular value decomposition and filtering
+ ! and matrices inversion to compute coefficient for the monomial
+ ! closure
USE prec_const
USE parallel
implicit none
-#ifdef TEST_SVD
! LOCAL VARIABLES
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: A_buff,Bf,Br ! buffer and full/reduced rebuilt matrices
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: U ! basis
@@ -11,19 +14,116 @@ module DLRA
INTEGER :: lda, ldu, ldvt, info, lwork, i, m,n, nsv_filter
COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: work
REAL(xp), DIMENSION(:), ALLOCATABLE :: rwork
- PUBLIC :: init_DLRA, filter_sv_moments_ky_pj, test_SVD
+
+ REAL(xp), DIMENSION(:), ALLOCATABLE :: ln ! Laguerre coefficients
+ REAL(xp), DIMENSION(:), ALLOCATABLE :: hn ! Hermite coefficients
+#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
CONTAINS
+ !--------------Routines to compute coeff for monomial truncation----------
+ ! Invert an upper triangular matrix
+ SUBROUTINE invert_utr_matrix(U,invU)
+ IMPLICIT NONE
+ REAL(xp), DIMENSION(:,:), INTENT(IN) :: U ! Matrix to invert
+ REAL(xp), DIMENSION(:,:), INTENT(OUT) :: invU! Result
+ END SUBROUTINE
+
+ SUBROUTINE build_hermite_matrix(d,U)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: d
+ REAL(xp), DIMENSION(d,d), INTENT(OUT) :: U
+ END SUBROUTINE
- SUBROUTINE init_DLRA(m_,n_)
+ SUBROUTINE build_laguerre_matrix(d,U)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: d
+ REAL(xp), DIMENSION(d,d), INTENT(OUT) :: U
+
+
+ END SUBROUTINE
+
+ ! 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)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: n ! Polyn. Deg. and index
+ INTEGER, DIMENSION(n+1) :: lnm2, lnm1
+ INTEGER :: k, e
+ ALLOCATE(ln(n+1))
+ IF (n == 0) THEN
+ ln = 1
+ ELSEIF (n == 1) THEN
+ ln = [-1, 1]
+ ELSE
+ lnm2 = 0.0
+ lnm2(n+1) = 1
+ lnm1 = 0.0
+ lnm1(n) = -1
+ lnm1(n+1) = 1
+ DO k = 2, n
+ ln = 0.0
+ DO e = n-k+1, n
+ ln(e) = (2*k-1)*lnm1(e) - lnm1(e+1) + (1-k)*lnm2(e)
+ END DO
+ ln(n+1) = (2*k-1)*lnm1(n+1) + (1-k)*lnm2(n+1)
+ ln = ln/k
+ 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
+
+ !--------------Routines to test singular value filtering -----------------
+#ifdef TEST_SVD
+ SUBROUTINE init_CLA(m_,n_)
USE basic
IMPLICIT NONE
! ARGUMENTS
INTEGER, INTENT(IN) :: m_,n_ ! dimensions of the input array
! read the input
INTEGER :: lun = 90 ! File duplicated from STDIN
- NAMELIST /DLRA/ nsv_filter
- READ(lun,dlra)
+ NAMELIST /CLA/ nsv_filter
+ READ(lun,CLA)
m = m_
n = n_
info = 1
@@ -46,14 +146,14 @@ CONTAINS
lwork = 2*CEILING(REAL(work(1)))
DEALLOCATE(work)
ALLOCATE(work(lwork))
- END SUBROUTINE init_DLRA
+ END SUBROUTINE init_CLA
SUBROUTINE filter_sv_moments_ky_pj
USE fields, ONLY: moments
USE grid, ONLY: total_nky, total_np, total_nj, ngp,ngj,ngz,&
local_np, local_nj, local_nz, local_nkx, local_nky, local_na
USE time_integration, ONLY: updatetlevel
- USE basic, ONLY: start_chrono, stop_chrono, chrono_DLRA
+ USE basic, ONLY: start_chrono, stop_chrono, chrono_CLA
IMPLICIT NONE
! Arguments
@@ -63,7 +163,7 @@ CONTAINS
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: moments_gky_lpj ! global ky, local pj data
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: moments_gky_gpj ! full gathered data for SVD (input of SVD)
INTEGER :: ia,ix,iz, m,n, ip, ij, iy
- CALL start_chrono(chrono_DLRA)
+ CALL start_chrono(chrono_CLA)
ALLOCATE(moments_lky_lpj(local_nky,local_np*local_nj))
ALLOCATE(moments_gky_lpj(total_nky,local_np*local_nj))
ALLOCATE(moments_gky_gpj(total_nky,total_np*total_nj))
@@ -121,7 +221,7 @@ CONTAINS
ENDDO
ENDDO
ENDDO
- CALL stop_chrono(chrono_DLRA)
+ CALL stop_chrono(chrono_CLA)
END SUBROUTINE filter_sv_moments_ky_pj
SUBROUTINE filter_singular_value(A)
@@ -142,7 +242,7 @@ CONTAINS
DO i=1,MIN(m,n)-nsv_filter
Sr(i) = Sf(i)
ENDDO
- ELSE ! do not filter if nsv_filter<0
+ ELSE ! DO not filter IF nsv_filter<0
Sr = Sf
ENDIF
! Reconstruct A from its reduced SVD
@@ -167,7 +267,7 @@ CONTAINS
END FUNCTION diagmat
SUBROUTINE test_svd
- ! Specify the dimensions of the input matrix A
+ ! SpecIFy the dimensions of the input matrix A
INTEGER :: m,n
! Declare the input matrix A and reconstructed matrix B
@@ -240,4 +340,4 @@ CONTAINS
stop
END SUBROUTINE test_svd
#endif
-end module DLRA
+END module CLA
diff --git a/src/auxval.F90 b/src/auxval.F90
index 1a043a31..5e69b5eb 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -13,7 +13,7 @@ subroutine auxval
USE parallel, ONLY: init_parallel_var, my_id, num_procs, num_procs_p, num_procs_z, num_procs_ky, rank_p, rank_ky, rank_z
USE processing, ONLY: init_process
#ifdef TEST_SVD
- USE DLRA, ONLY: init_DLRA
+ USE CLA, ONLY: init_CLA
#endif
IMPLICIT NONE
INTEGER :: i_, ierr
@@ -44,7 +44,7 @@ subroutine auxval
#ifdef TEST_SVD
! If we want to test SVD decomposition etc.
- CALL init_DLRA(local_nky,local_np*local_nj)
+ CALL init_CLA(local_nky,local_np*local_nj)
#endif
!! Display parallel settings
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 7119944d..704aad2a 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -37,7 +37,7 @@ MODULE basic
chrono_diag, chrono_chck, chrono_step, chrono_clos, chrono_ghst, chrono_coll, chrono_napj, chrono_grad
#ifdef TEST_SVD
! A chrono for SVD tests
- type(chrono), PUBLIC, PROTECTED :: chrono_DLRA
+ type(chrono), PUBLIC, PROTECTED :: chrono_CLA
#endif
! This sets if the outputs is done through a large gather or using parallelization from futils
! it is recommended to set it to .true.
@@ -87,7 +87,7 @@ CONTAINS
chrono_diag%ttot = 0
chrono_step%ttot = 0
#ifdef TEST_SVD
- chrono_DLRA%ttot = 0
+ chrono_CLA%ttot = 0
#endif
END SUBROUTINE basic_data
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 32e2c137..a27cd910 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -36,7 +36,7 @@ SUBROUTINE set_closure_model
! set the evolve mom array
ALLOCATE(evolve_mom(local_np+ngp,local_nj+ngj))
SELECT CASE(hierarchy_closure)
- CASE('truncation')
+ CASE('truncation','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)) &
@@ -56,7 +56,7 @@ SUBROUTINE set_closure_model
! Set the nonlinear closure scheme (truncation of sum over n in Sapj)
ALLOCATE(nmaxarray(local_nj))
SELECT CASE(nonlinear_closure)
- CASE('truncation')
+ CASE('truncation','monomial')
IF(nmax .LT. 0) THEN
CALL speak("Set nonlinear truncation to anti Laguerre aliasing")
DO ij = 1,local_nj
@@ -75,6 +75,11 @@ SUBROUTINE set_closure_model
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
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index e1914a7a..f9e3aa8d 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -131,7 +131,7 @@ SUBROUTINE diagnose_full(kstep)
CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
#ifdef TEST_SVD
- CALL creatd(fidres, 0, (/0/), "/profiler/Tc_DLRA", "cumulative total DLRA computation time")
+ CALL creatd(fidres, 0, (/0/), "/profiler/Tc_CLA", "cumulative total CLA computation time")
#endif
! Grid info
CALL creatg(fidres, "/data/grid", "Grid data")
@@ -328,7 +328,7 @@ SUBROUTINE diagnose_0d
CALL append(fidres, "/profiler/Tc_nadiab", REAL(chrono_napj%ttot,dp),ionode=0)
CALL append(fidres, "/profiler/Tc_step", REAL(chrono_step%ttot,dp),ionode=0)
#ifdef TEST_SVD
- CALL append(fidres, "/profiler/Tc_DLRA", REAL(chrono_DLRA%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_CLA", REAL(chrono_CLA%ttot,dp),ionode=0)
#endif
CALL append(fidres, "/profiler/time", REAL(time,dp),ionode=0)
! Processing data
@@ -355,7 +355,7 @@ SUBROUTINE diagnose_2d
USE diagnostics_par
USE futils, ONLY: putarr, append
#ifdef TEST_SVD
- USE DLRA, ONLY: Sf
+ USE CLA, ONLY: Sf
#endif
IMPLICIT NONE
iframe2d=iframe2d+1
diff --git a/src/stepon.F90 b/src/stepon.F90
index 16a25bbd..af24f831 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -9,7 +9,7 @@ SUBROUTINE stepon
USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: xp
#ifdef TEST_SVD
- USE DLRA, ONLY: test_svd,filter_sv_moments_ky_pj
+ USE CLA, ONLY: test_svd,filter_sv_moments_ky_pj
#endif
IMPLICIT NONE
diff --git a/testcases/DLRA_zpinch/base_case/fort_00.90 b/testcases/DLRA_zpinch/base_case/fort_00.90
index 6caf3d73..07e2cd32 100644
--- a/testcases/DLRA_zpinch/base_case/fort_00.90
+++ b/testcases/DLRA_zpinch/base_case/fort_00.90
@@ -102,6 +102,6 @@
&TIME_INTEGRATION_PAR
numerical_scheme = 'RK4'
/
-&DLRA
+&CLA
nsv_filter = 0
/
diff --git a/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
index 1235695e..1fa68c8b 100644
--- a/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
+++ b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
@@ -102,6 +102,6 @@
&TIME_INTEGRATION_PAR
numerical_scheme = 'RK4'
/
-&DLRA
+&CLA
nsv_filter = 2
/
diff --git a/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
index a23d82d4..81ae89e5 100644
--- a/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
+++ b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
@@ -102,6 +102,6 @@
&TIME_INTEGRATION_PAR
numerical_scheme = 'RK4'
/
-&DLRA
+&CLA
nsv_filter = 6
/
--
GitLab