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

Testing of singular value filtering is operational

parent d157cc1c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 613 additions and 230 deletions
Makefile
+ 5
− 3
View file @ 48150a74
......@@ -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 $@
src/DLRA_mod.F90
+ 227
− 223
View file @ 48150a74
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
src/auxval.F90
+ 7
− 0
View file @ 48150a74
......@@ -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
......
src/basic_mod.F90
+ 7
− 0
View file @ 48150a74
......@@ -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
......
src/diagnose.F90
+ 37
− 1
View file @ 48150a74
......@@ -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
......
src/stepon.F90
+ 5
− 3
View file @ 48150a74
......@@ -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
......
testcases/DLRA_zpinch/base_case/fort_00.90 0 → 100644
+ 107
− 0
View file @ 48150a74
&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
/
testcases/DLRA_zpinch/base_case/gyacomo23_sp 0 → 120000
+ 1
− 0
View file @ 48150a74
../../../bin/gyacomo23_sp
\ No newline at end of file
testcases/DLRA_zpinch/base_case/gyacomo23_test_svd 0 → 120000
+ 1
− 0
View file @ 48150a74
../../../bin/gyacomo23_test_svd
\ No newline at end of file
testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 0 → 100644
+ 107
− 0
View file @ 48150a74
&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
/
testcases/DLRA_zpinch/nsv_filter_2/gyacomo23_test_svd 0 → 120000
+ 1
− 0
View file @ 48150a74
../../../bin/gyacomo23_test_svd
\ No newline at end of file
testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 0 → 100644
+ 107
− 0
View file @ 48150a74
&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
/
testcases/DLRA_zpinch/nsv_filter_6/gyacomo23_test_svd 0 → 120000
+ 1
− 0
View file @ 48150a74
../../../bin/gyacomo23_test_svd
\ No newline at end of file
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