From 48150a743fecb00e634c8980ba52260a5e4fdd79 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 6 Apr 2023 15:52:17 +0200
Subject: [PATCH] Testing of singular value filtering is operational
---
Makefile | 8 +-
src/DLRA_mod.F90 | 450 +++++++++---------
src/auxval.F90 | 7 +
src/basic_mod.F90 | 7 +
src/diagnose.F90 | 38 +-
src/stepon.F90 | 8 +-
testcases/DLRA_zpinch/base_case/fort_00.90 | 107 +++++
testcases/DLRA_zpinch/base_case/gyacomo23_sp | 1 +
.../DLRA_zpinch/base_case/gyacomo23_test_svd | 1 +
testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 | 107 +++++
.../nsv_filter_2/gyacomo23_test_svd | 1 +
testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 | 107 +++++
.../nsv_filter_6/gyacomo23_test_svd | 1 +
13 files changed, 613 insertions(+), 230 deletions(-)
create mode 100644 testcases/DLRA_zpinch/base_case/fort_00.90
create mode 120000 testcases/DLRA_zpinch/base_case/gyacomo23_sp
create mode 120000 testcases/DLRA_zpinch/base_case/gyacomo23_test_svd
create mode 100644 testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
create mode 120000 testcases/DLRA_zpinch/nsv_filter_2/gyacomo23_test_svd
create mode 100644 testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
create mode 120000 testcases/DLRA_zpinch/nsv_filter_6/gyacomo23_test_svd
diff --git a/Makefile b/Makefile
index a0e17d3b..8e713893 100644
--- a/Makefile
+++ b/Makefile
@@ -61,7 +61,8 @@ gdebug: compile
# test SVD
test_svd: F90 = mpif90
-test_svd: F90FLAGS = -DTEST_SVD -g -traceback -ftrapuv -warn all -debug all
+# test_svd: F90FLAGS = -DTEST_SVD -g -traceback -ftrapuv -warn all -debug all
+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
@@ -123,7 +124,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)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/DLRA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
$(OBJDIR)/basic_mod.o : src/basic_mod.F90 \
@@ -309,5 +310,6 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
$(OBJDIR)/DLRA_mod.o : src/DLRA_mod.F90 \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(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 $@
diff --git a/src/DLRA_mod.F90 b/src/DLRA_mod.F90
index c6a39166..3be00e96 100644
--- a/src/DLRA_mod.F90
+++ b/src/DLRA_mod.F90
@@ -1,239 +1,243 @@
module DLRA
- USE prec_const
- implicit none
-
- PUBLIC :: filter_singular_value_ky_pj, test_SVD
-
- CONTAINS
-
- SUBROUTINE filter_singular_value_ky_pj(nsv,array_ky_pj)
- IMPLICIT NONE
- ! ARGUMENTS
- INTEGER, INTENT(IN) :: nsv ! number of singular values to keep
- COMPLEX(xp), DIMENSION(:,:), INTENT(INOUT) :: array_ky_pj ! Array to filter
- !
-
- ! Singular value decomposition
- ! CALL SVD(array_ky_pj,singular_values)
- END SUBROUTINE
-
- SUBROUTINE test_svd
+ USE prec_const
+ USE parallel
+ implicit none
#ifdef TEST_SVD
- ! Program to perform Singular Value Decomposition (SVD)
- ! using LAPACK library
-
- ! Specify the dimensions of the input matrix A
- INTEGER, PARAMETER :: m = 3, n = 2
-
- ! Declare the input matrix A
- COMPLEX, DIMENSION(m,n) :: A
-
- ! OUTPUT
- COMPLEX, DIMENSION(m,m) :: U
- REAL, DIMENSION(MIN(m,n)) :: S
- COMPLEX, DIMENSION(n,n) :: VT
-
- ! local variables
- INTEGER :: lda, ldu, ldvt, info, lwork, i, j
- COMPLEX, DIMENSION(:), ALLOCATABLE :: work
- REAL, DIMENSION(:), ALLOCATABLE :: rwork
-
- ! Set the leading dimensions for the input and output arrays
- lda = MAX(1, m)
- ldu = MAX(1, m)
- ldvt = MAX(1, n)
-
- ! Define the input matrix A
- A = RESHAPE((/ (1.0,0.1), (2.0,0.2), (3.0,0.3), (4.0,0.4), (5.0,0.5), (6.0,0.6) /), SHAPE(A))
-
- ! Print the input matrix A
- WRITE(*,*) 'Input matrix A = '
- DO i = 1, m
- WRITE(*,*) ('(',REAL(A(i,j)), AIMAG(A(i,j)),')', j=1,n)
- END DO
-
- ALLOCATE(work(5*n), rwork(5*n))
- ! Compute the optimal workspace size
- lwork = -1
- CALL CGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
-
- ! Allocate memory for the workspace arrays
- lwork = CEILING(REAL(work(1)), KIND=SELECTED_REAL_KIND(1, 6))
-
- ! Compute the SVD of A using the LAPACK subroutine CGESVD
- CALL CGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
-
- ! Print the results
- WRITE(*,*) 'U = '
- DO i = 1, m
- WRITE(*,*) ('(',REAL(U(i,j)), AIMAG(U(i,j)),')', j=1,m)
- END DO
- WRITE(*,*)
- WRITE(*,*) 'S = '
- WRITE(*,'(2F8.3)') (S(i), i=1,n)
- WRITE(*,*)
- WRITE(*,*) 'VT = '
- DO i = 1, n
- WRITE(*,*) ('(',REAL(VT(i,j)), AIMAG(VT(i,j)),')', j=1,n)
- END DO
-
- ! Reconstruct A from its SVD
- A = MATMUL(U, MATMUL(diagmat(S,m,n), VT))
-
- ! Print the reconstructed matrix A
- WRITE(*,*) 'Reconstructed matrix A = '
- DO i = 1, m
- WRITE(*,*) ('(',REAL(A(i,j)), AIMAG(A(i,j)),')', j=1,n)
- END DO
- stop
+ ! LOCAL VARIABLES
+ 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
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: VT ! basis
+ 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
+
+CONTAINS
+
+ SUBROUTINE init_DLRA(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)
+ m = m_
+ n = n_
+ info = 1
+ ! Allocate the matrices
+ ALLOCATE(A_buff(m,n),Bf(m,n),Br(m,n))
+ ALLOCATE(U(m,m),Sf(MIN(m,n)),Sr(MIN(m,n)),VT(n,n))
+ ! Set the leading dimensions for the input and output arrays
+ lda = MAX(1, m)
+ ldu = MAX(1, m)
+ ldvt = MAX(1, n)
+ ALLOCATE(work(5*n), rwork(5*MIN(m,n)))
+ ! Compute the optimal workspace size
+ lwork = -1
+#ifdef SINGLE_PRECISION
+ CALL CGESVD('A', 'A', m, n, A_buff, lda, Sf, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#else
+ CALL ZGESVD('A', 'A', m, n, A_buff, lda, Sf, U, ldu, VT, ldvt, work, lwork, rwork, info)
#endif
- END SUBROUTINE test_svd
-
- ! SUBROUTINE test_svd
- ! ! Program to perform Singular Value Decomposition (SVD)
- ! ! using LAPACK library
-
- ! ! Specify the dimensions of the input matrix A
- ! INTEGER, PARAMETER :: m = 3, n = 2
- ! INTEGER, PARAMETER :: lda = m
-
- ! ! Declare the input matrix A
- ! REAL, DIMENSION(lda,n) :: A
-
- ! ! Specify the dimensions of the output matrices
- ! INTEGER, PARAMETER :: ldu = m, ldvt = n
- ! INTEGER, PARAMETER :: lwork = 5*n
- ! REAL, DIMENSION(ldu,m) :: U
- ! REAL, DIMENSION(n) :: S
- ! REAL, DIMENSION(ldvt,n) :: VT
- ! REAL, DIMENSION(lwork) :: work
- ! INTEGER :: info,i,j
-
- ! ! Define the input matrix A
- ! A = RESHAPE((/ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 /), SHAPE(A))
-
- ! ! Compute the SVD of A using the LAPACK subroutine SGESVD
- ! CALL SGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, info)
-
- ! ! Print the results
- ! WRITE(*,*) 'U = '
- ! DO i = 1, m
- ! WRITE(*,'(6F8.3)') (U(i,j), j=1,m)
- ! END DO
- ! WRITE(*,*)
- ! WRITE(*,*) 'S = '
- ! WRITE(*,'(2F8.3)') (S(i), i=1,n)
- ! WRITE(*,*)
- ! WRITE(*,*) 'VT = '
- ! DO i = 1, n
- ! WRITE(*,'(6F8.3)') (VT(i,j), j=1,n)
- ! END DO
-
- ! ! Reconstruct A from its SVD
- ! A = MATMUL(U, MATMUL(diagmat(S,m,n), TRANSPOSE(VT)))
-
- ! ! Print the reconstructed matrix A
- ! WRITE(*,*) 'Reconstructed matrix A = '
- ! DO i = 1, m
- ! WRITE(*,'(2X, 3F8.3)') (A(i,j), j=1,n)
- ! END DO
- ! stop
- ! END SUBROUTINE test_svd
-
- ! SUBROUTINE test_svd
- ! IMPLICIT NONE
- ! INTEGER, PARAMETER :: m = 3, n = 2
- ! COMPLEX(xp), DIMENSION(m, n) :: A
- ! COMPLEX(xp), DIMENSION(m, m) :: U
- ! COMPLEX(xp), DIMENSION(n, n) :: VT
- ! REAL(xp), DIMENSION(MIN(m,n)) :: S
- ! INTEGER :: i, j
-
- ! ! Initialize A
- ! A = RESHAPE((/ (1.0_xp, 2.0_xp), (3.0_xp, 4.0_xp), (5.0_xp, 6.0_xp),&
- ! (1.5_xp, 2.5_xp), (3.5_xp, 4.5_xp), (5.5_xp, 6.5_xp) /), [m, n])
-
- ! ! Print input
- ! WRITE(*,*) "A = "
- ! DO i = 1, m
- ! WRITE(*,"(3F8.3)") (REAL(A(i,j)), AIMAG(A(i,j)), j=1,n)
- ! END DO
- ! WRITE(*,*)
-
- ! ! Call the SVD subroutine
- ! CALL svd(A, U, S, VT, m, n)
-
- ! ! Print the resultas
- ! WRITE(*,*) "U = "
- ! DO i = 1, m
- ! WRITE(*,"(3F8.3)") (REAL(U(i,j)), AIMAG(U(i,j)), j=1,m)
- ! END DO
- ! WRITE(*,*)
-
- ! WRITE(*,*) "S = ", S
- ! WRITE(*,*)
-
- ! WRITE(*,*) "VT = "
- ! DO i = 1, n
- ! WRITE(*,"(3F8.3)") (REAL(VT(i,j)), AIMAG(VT(i,j)), j=1,n)
- ! END DO
- ! END SUBROUTINE test_svd
-
-
- SUBROUTINE svd(A, U, S, VT, m, n)
+ ! Allocate memory for the workspace arrays
+ lwork = 2*CEILING(REAL(work(1)))
+ DEALLOCATE(work)
+ ALLOCATE(work(lwork))
+ END SUBROUTINE init_DLRA
+
+ 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
IMPLICIT NONE
- ! INPUT
- COMPLEX(xp), DIMENSION(m,n), INTENT(IN) :: A
- INTEGER, INTENT(IN) :: m, n
- ! OUTPUT
- COMPLEX(xp), DIMENSION(m,m), INTENT(OUT) :: U
- REAL(xp), DIMENSION(MIN(m,n)), INTENT(OUT) :: S
- COMPLEX(xp), DIMENSION(n,n), INTENT(OUT) :: VT
- ! local variables
- INTEGER :: lda, ldu, ldvt, info, lwork
- COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: work
- REAL(xp), DIMENSION(:), ALLOCATABLE :: rwork
-#ifdef LAPACKDIR
- ! Set the leading dimensions for the input and output arrays
- lda = MAX(1, m)
- ldu = MAX(1, m)
- ldvt = MAX(1, n)
-
- ! Compute the optimal workspace size
- lwork = -1
- CALL ZGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
-
- ! Allocate memory for the workspace arrays
- lwork = CEILING(REAL(work(1)), KIND=SELECTED_REAL_KIND(1, 6))
- ALLOCATE(work(lwork), rwork(5*MIN(m,n)))
-
- ! Compute the SVD
- CALL ZGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
-
- ! Free the workspace arrays
- DEALLOCATE(work, rwork)
+
+ ! Arguments
+ INTEGER :: nsv_filter ! number of singular values to keep
+ ! Local variables
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: moments_lky_lpj ! local ky and local pj data
+ 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)
+ 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))
+ DO iz = 1+ngz/2,local_nz+ngz/2
+ DO ix = 1,local_nkx
+ DO ia = 1,local_na
+ ! Build the local slice explicitely
+ DO ip = 1,local_np
+ DO ij = 1,local_nj
+ DO iy = 1,local_nky
+ moments_lky_lpj(iy,local_np*(ij-1)+ip) = moments(ia,ip+ngp/2,ij+ngj/2,iy,ix,iz,updatetlevel)
+ ENDDO
+ ENDDO
+ ENDDO
+ ! Gather ky data
+ IF(num_procs_ky .GT. 1) THEN
+ ! MPI communication
+ ELSE
+ moments_gky_lpj = moments_lky_lpj
+ ENDIF
+ ! Gather p data
+ IF(num_procs_p .GT. 1) THEN
+ ! MPI communication
+ ELSE
+ moments_gky_gpj = moments_gky_lpj
+ ENDIF
+ ! The process 0 performs the SVD
+ IF(my_id .EQ. 0) THEN
+ m = total_nky
+ n = total_np*total_nj
+ nsv_filter = -1
+ CALL filter_singular_value(moments_gky_gpj)
+ ENDIF
+ ! Distribute ky data
+ IF(num_procs_ky .GT. 1) THEN
+ ! MPI communication
+ ELSE
+ moments_gky_lpj = moments_gky_gpj
+ ENDIF
+ ! Distribute p data
+ IF(num_procs_p .GT. 1) THEN
+ ! MPI communication
+ ELSE
+ moments_lky_lpj = moments_gky_lpj
+ ENDIF
+ ! Put back the data into the moments array
+ DO ip = 1,local_np
+ DO ij = 1,local_nj
+ DO iy = 1,local_nky
+ moments(ia,ip+ngp/2,ij+ngj/2,iy,ix,iz,updatetlevel) = moments_lky_lpj(iy,local_np*(ij-1)+ip)
+ ENDDO
+ ENDDO
+ ENDDO
+
+ ENDDO
+ ENDDO
+ ENDDO
+ CALL stop_chrono(chrono_DLRA)
+ END SUBROUTINE filter_sv_moments_ky_pj
+
+ SUBROUTINE filter_singular_value(A)
+ IMPLICIT NONE
+ ! ARGUMENTS
+ COMPLEX(xp), DIMENSION(:,:), INTENT(INOUT) :: A ! Array to filter
+ ! copy of the input since it is changed in the svd procedure
+ A_buff = A;
+#ifdef SINGLE_PRECISION
+ CALL CGESVD('A', 'A', m, n, A_buff, lda, Sf, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#else
+ CALL ZGESVD('A', 'A', m, n, A_buff, lda, Sf, U, ldu, VT, ldvt, work, lwork, rwork, info)
#endif
- END SUBROUTINE svd
-
-
- FUNCTION diagmat(S, m, n) RESULT(D)
+ IF(info.GT.0) print*,"SVD did not converge, info=",info
+ ! Filter the values
+ Sr = 0._xp
+ IF (nsv_filter .GT. 0) THEN
+ DO i=1,MIN(m,n)-nsv_filter
+ Sr(i) = Sf(i)
+ ENDDO
+ ELSE ! do not filter if nsv_filter<0
+ Sr = Sf
+ ENDIF
+ ! Reconstruct A from its reduced SVD
+ A = MATMUL(U, MATMUL(diagmat(Sr,m,n), VT)) ! reduced
+ END SUBROUTINE filter_singular_value
+
+ FUNCTION diagmat(S, m, n) RESULT(D)
IMPLICIT NONE
! INPUT
- REAL, DIMENSION(:), INTENT(IN) :: S
+ REAL(xp), DIMENSION(:), INTENT(IN) :: S
INTEGER, INTENT(IN) :: m, n
! OUTPUT
- COMPLEX, DIMENSION(m,n) :: D
+ COMPLEX(xp), DIMENSION(m,n) :: D
! Local variables
- INTEGER :: i, j
-
+ INTEGER :: i
! Initialize the output array to zero
D = 0.0
-
! Fill the diagonal elements of the output array with the input vector S
DO i = 1, MIN(m,n)
- D(i,i) = S(i)
+ D(i,i) = S(i)
END DO
-
- END FUNCTION diagmat
+ END FUNCTION diagmat
-end module DLRA
\ No newline at end of file
+ SUBROUTINE test_svd
+ ! Specify the dimensions of the input matrix A
+ INTEGER :: m,n
+
+ ! Declare the input matrix A and reconstructed matrix B
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: A,B,A_buff
+
+ ! OUTPUT
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: U
+ REAL(xp), DIMENSION(:), ALLOCATABLE :: S
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: VT
+ ! local variables
+ INTEGER :: lda, ldu, ldvt, info, lwork, i, j
+ COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: work
+ REAL(xp), DIMENSION(:), ALLOCATABLE :: rwork
+
+ m = 3
+ n = 2
+ ALLOCATE(A_buff(m,n),A(m,n),B(m,n),U(m,m),S(MIN(m,n)),VT(n,n))
+ ! Set the leading dimensions for the input and output arrays
+ lda = MAX(1, m)
+ ldu = MAX(1, m)
+ ldvt = MAX(1, n)
+
+ ! Define the input matrix A
+ A = RESHAPE((/ (1._xp,0.1_xp), (2._xp,0.2_xp), (3._xp,0.3_xp), (4._xp,0.4_xp), (5._xp,0.5_xp), (6._xp,0.6_xp) /), SHAPE(A))
+ ! copy of the input since it is changed in the svd procedure
+ A_buff = A;
+ ! Print the input matrix A
+ WRITE(*,*) 'Input matrix A = '
+ DO i = 1, m
+ WRITE(*,*) ('(',REAL(A(i,j)), AIMAG(A(i,j)),')', j=1,n)
+ END DO
+
+ ! CALL svd(A, U, S, VT, m, n)
+ ALLOCATE(work(5*n), rwork(5*n))
+ ! Compute the optimal workspace size
+ lwork = -1
+#ifdef SINGLE_PRECISION
+ CALL CGESVD('A', 'A', m, n, A_buff, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#else
+ CALL ZGESVD('A', 'A', m, n, A_buff, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#endif
+
+ ! Allocate memory for the workspace arrays
+ lwork = CEILING(REAL(work(1)), KIND=SELECTED_REAL_KIND(1, 6))
+
+ ! Compute the SVD of A using the LAPACK subroutine CGESVD
+#ifdef SINGLE_PRECISION
+ CALL CGESVD('A', 'A', m, n, A_buff, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#else
+ CALL ZGESVD('A', 'A', m, n, A_buff, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
+#endif
+ ! Print the results
+ ! WRITE(*,*) 'U = '
+ ! DO i = 1, m
+ ! WRITE(*,*) ('(',REAL(U(i,j)), AIMAG(U(i,j)),')', j=1,m)
+ ! END DO
+ ! WRITE(*,*)
+ WRITE(*,*) 'S = '
+ WRITE(*,'(2F8.3)') (S(i), i=1,n)
+ WRITE(*,*)
+ ! WRITE(*,*) 'VT = '
+ ! DO i = 1, n
+ ! WRITE(*,*) ('(',REAL(VT(i,j)), AIMAG(VT(i,j)),')', j=1,n)
+ ! END DO
+
+ ! Reconstruct A from its SVD
+ B = MATMUL(U, MATMUL(diagmat(S,m,n), VT))
+ ! Print the error with the intput matrix
+ WRITE(*,*) '||A-USVT||=', sum(abs(A-B))
+ stop
+ END SUBROUTINE test_svd
+#endif
+end module DLRA
diff --git a/src/auxval.F90 b/src/auxval.F90
index fb8e13ca..136e4bc0 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -12,6 +12,9 @@ subroutine auxval
USE closure, ONLY: set_closure_model, hierarchy_closure
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
+#endif
IMPLICIT NONE
INTEGER :: i_, ierr
@@ -42,6 +45,10 @@ subroutine auxval
CALL set_closure_model ! set the closure scheme in use
+#ifdef TEST_SVD
+ CALL init_DLRA(local_nky,local_np*local_nj)
+#endif
+
!! Display parallel settings
CALL mpi_barrier(MPI_COMM_WORLD, ierr)
DO i_ = 0,num_procs-1
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 0784de19..2aa6dfd8 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -40,6 +40,10 @@ MODULE basic
type(chrono), PUBLIC, PROTECTED :: chrono_runt, chrono_mrhs, chrono_advf, chrono_pois, chrono_sapj,&
chrono_diag, chrono_chck, chrono_step, chrono_clos, chrono_ghst, chrono_coll, chrono_napj, chrono_grad
+#ifdef TEST_SVD
+ type(chrono), PUBLIC, PROTECTED :: chrono_DLRA
+#endif
+
LOGICAL, PUBLIC, PROTECTED :: GATHERV_OUTPUT = .true.
PUBLIC :: allocate_array, basic_outputinputs,basic_data,&
@@ -85,6 +89,9 @@ CONTAINS
chrono_chck%ttot = 0
chrono_diag%ttot = 0
chrono_step%ttot = 0
+#ifdef TEST_SVD
+ chrono_DLRA%ttot = 0
+#endif
END SUBROUTINE basic_data
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index c3be6863..e9f89c58 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -127,6 +127,9 @@ SUBROUTINE diagnose_full(kstep)
CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time")
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")
+#endif
! Grid info
CALL creatg(fidres, "/data/grid", "Grid data")
CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0)
@@ -169,6 +172,15 @@ SUBROUTINE diagnose_full(kstep)
ENDIF
CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
END IF
+ ! var2d group (gyro transport)
+ IF (nsave_0d .GT. 0) THEN
+ CALL creatg(fidres, "/data/var2d", "2d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
+#ifdef TEST_SVD
+ CALL creatg(fidres, "/data/var2d/sv_ky_pj", "singular values of the moment ky/pj cut")
+#endif
+ ENDIF
! var3d group (phi,psi, fluid moments, Ni00, Napjz)
IF (nsave_3d .GT. 0) THEN
CALL creatg(fidres, "/data/var3d", "3d profiles")
@@ -221,6 +233,12 @@ SUBROUTINE diagnose_full(kstep)
CALL diagnose_0d
END IF
END IF
+ ! 2.2 2d profiles
+ IF (nsave_2d .GT. 0) THEN
+ IF (MOD(cstep, nsave_2d) == 0) THEN
+ CALL diagnose_2d
+ ENDIF
+ ENDIF
! 2.3 3d profiles
IF (nsave_3d .GT. 0) THEN
IF (MOD(cstep, nsave_3d) == 0) THEN
@@ -274,7 +292,6 @@ 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 creatd(fidres, 0, (/0/), "/profiler/Tc_DLRA", "cumulative total DLRA computation time")
CALL append(fidres, "/profiler/Tc_DLRA", REAL(chrono_DLRA%ttot,dp),ionode=0)
#endif
CALL append(fidres, "/profiler/time", REAL(time,dp),ionode=0)
@@ -296,6 +313,25 @@ SUBROUTINE diagnose_0d
ENDIF
END SUBROUTINE diagnose_0d
+SUBROUTINE diagnose_2d
+ USE prec_const
+ USE basic
+ USE diagnostics_par
+ USE futils, ONLY: putarr, append
+#ifdef TEST_SVD
+ USE DLRA, ONLY: Sf
+#endif
+ IMPLICIT NONE
+ CHARACTER(50) :: dset_name
+ iframe2d=iframe2d+1
+ CALL append(fidres,"/data/var2d/time", REAL(time,dp), ionode=0)
+ CALL append(fidres,"/data/var2d/cstep",REAL(cstep,dp),ionode=0)
+#ifdef TEST_SVD
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var2d/sv_ky_pj/", iframe2d
+ CALL putarr(fidres, dset_name, Sf, ionode=0)
+#endif
+END SUBROUTINE
+
SUBROUTINE diagnose_3d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd, putarr
diff --git a/src/stepon.F90 b/src/stepon.F90
index 71ad932a..930a6b2c 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -8,7 +8,9 @@ SUBROUTINE stepon
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: xp
- USE DLRA, ONLY: test_svd
+#ifdef TEST_SVD
+ USE DLRA, ONLY: test_svd,filter_sv_moments_ky_pj
+#endif
IMPLICIT NONE
INTEGER :: num_step, ierr
@@ -56,9 +58,9 @@ SUBROUTINE stepon
CALL stop_chrono(chrono_chck)
!! TEST SINGULAR VALUE DECOMPOSITION
- ! CALL filter_singular_value_ky_pj(nsv,moments)
#ifdef TEST_SVD
- CALL test_svd
+ ! CALL test_svd
+ CALL filter_sv_moments_ky_pj
#endif
IF( nlend ) EXIT ! exit do loop
diff --git a/testcases/DLRA_zpinch/base_case/fort_00.90 b/testcases/DLRA_zpinch/base_case/fort_00.90
new file mode 100644
index 00000000..6caf3d73
--- /dev/null
+++ b/testcases/DLRA_zpinch/base_case/fort_00.90
@@ -0,0 +1,107 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 500
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 4
+ jmax = 2
+ Nx = 64
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.1
+ dtsave_1d = -1
+ dtsave_2d = 0.1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ !hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
+&DLRA
+ nsv_filter = 0
+/
diff --git a/testcases/DLRA_zpinch/base_case/gyacomo23_sp b/testcases/DLRA_zpinch/base_case/gyacomo23_sp
new file mode 120000
index 00000000..db5b14c7
--- /dev/null
+++ b/testcases/DLRA_zpinch/base_case/gyacomo23_sp
@@ -0,0 +1 @@
+../../../bin/gyacomo23_sp
\ No newline at end of file
diff --git a/testcases/DLRA_zpinch/base_case/gyacomo23_test_svd b/testcases/DLRA_zpinch/base_case/gyacomo23_test_svd
new file mode 120000
index 00000000..d0950bb2
--- /dev/null
+++ b/testcases/DLRA_zpinch/base_case/gyacomo23_test_svd
@@ -0,0 +1 @@
+../../../bin/gyacomo23_test_svd
\ No newline at end of file
diff --git a/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
new file mode 100644
index 00000000..1235695e
--- /dev/null
+++ b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
@@ -0,0 +1,107 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 500
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 4
+ jmax = 2
+ Nx = 64
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.1
+ dtsave_1d = -1
+ dtsave_2d = 0.1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ !hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
+&DLRA
+ nsv_filter = 2
+/
diff --git a/testcases/DLRA_zpinch/nsv_filter_2/gyacomo23_test_svd b/testcases/DLRA_zpinch/nsv_filter_2/gyacomo23_test_svd
new file mode 120000
index 00000000..d0950bb2
--- /dev/null
+++ b/testcases/DLRA_zpinch/nsv_filter_2/gyacomo23_test_svd
@@ -0,0 +1 @@
+../../../bin/gyacomo23_test_svd
\ No newline at end of file
diff --git a/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
new file mode 100644
index 00000000..a23d82d4
--- /dev/null
+++ b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
@@ -0,0 +1,107 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 500
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 4
+ jmax = 2
+ Nx = 64
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.1
+ dtsave_1d = -1
+ dtsave_2d = 0.1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ !hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
+&DLRA
+ nsv_filter = 6
+/
diff --git a/testcases/DLRA_zpinch/nsv_filter_6/gyacomo23_test_svd b/testcases/DLRA_zpinch/nsv_filter_6/gyacomo23_test_svd
new file mode 120000
index 00000000..d0950bb2
--- /dev/null
+++ b/testcases/DLRA_zpinch/nsv_filter_6/gyacomo23_test_svd
@@ -0,0 +1 @@
+../../../bin/gyacomo23_test_svd
\ No newline at end of file
--
GitLab