diff --git a/Makefile b/Makefile index 8c27c4c..b107bd1 100755 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ FFILES = rhs.f90 vis.f90 fluid_time_step.f90 init_fields.f90 \ sponge.f90 fft_unit_test.f90 draw_plate.f90 draw_sphere.f90 \ rotation_matrices.f90 add_channel.f90 add_cavity.f90 init_scalar.f90 dry_run.f90 \ noncircular_cylinder.f90 draw_flexible_plate.f90 \ - runtime_backuping.f90 io_test.f90 POD.f90 + runtime_backuping.f90 io_test.f90 POD.f90 post_force.f90 ifndef NOHDF5 # Case WITH HDF5 (all machines except earth simulator) diff --git a/src/ini_files_parser.f90 b/src/ini_files_parser.f90 index fcbdd29..31f719d 100644 --- a/src/ini_files_parser.f90 +++ b/src/ini_files_parser.f90 @@ -59,46 +59,65 @@ module module_ini_files_parser ! note: array is assumed-shape and its size defines what we try to read ! --> MPI wrapper in the MPI parser module !----------------------------------------------------------------------------- - subroutine read_array_from_ascii_file(file, array, n_header) - implicit none - character(len=*), intent(in) :: file - integer, intent(in) :: n_header - real(kind=pr), intent(inout) :: array (1:,1:) - integer :: nlines, ncols, i, io_error - character(len=maxcolumns) :: dummy - character(len=16) :: fmt - character(len=3) :: ncols_str - - ! check if the specified file exists - call check_file_exists( file ) - - nlines = size(array,1) - ncols = size(array,2) - - write(*,'(80("-"))') - write(*,'("INFO: reading ",i8," lines with ",i8," colums from ",A)') nlines, ncols, file - - ! set up format string - write(ncols_str,'(i3.3)') ncols - fmt = '('//ncols_str//'(es12.4,1x))' - - io_error = 0 - i = 0 - - open(unit=14,file=trim(adjustl(file)),action='read',status='old') - do while (io_error==0) - ! read a line from file - read (14,'(A)',iostat=io_error) dummy - i = i + 1 - ! if we're past the header AND the read worked (i.e. not end of file) - if (i > n_header .and. io_error==0) then - read(dummy,*) array(i-n_header,:) + subroutine read_array_from_ascii_file(file, array, n_header, user_separator) + implicit none + character(len=*), intent(in) :: file + integer, intent(in) :: n_header + real(kind=pr), intent(inout) :: array (1:,1:) + character(len=1), optional, intent(in) :: user_separator + integer :: nlines, ncols, i, io_error, k + character(len=maxcolumns) :: dummy + character(len=16) :: fmt + character(len=3) :: ncols_str + character(len=1) :: sep + + ! check if the specified file exists + call check_file_exists( file ) + + nlines = size(array,1) + ncols = size(array,2) + + if (present(user_separator)) then + sep = user_separator + else + sep = " " endif - enddo - close (14) - write(*,'("Done reading.")') - write(*,'(80("-"))') + write(*,'(80("-"))') + write(*,'("INFO: reading ",i8," lines with ",i8," colums from ",A)') nlines, ncols, file + write(*,'("INFO: separator is -->",A1,"<--")') sep + + ! set up format string + write(ncols_str,'(i3.3)') ncols + fmt = '('//ncols_str//'(es12.4,1x))' + + io_error = 0 + i = 0 + + open(unit=14,file=trim(adjustl(file)),action='read',status='old') + do while (io_error==0) + ! read a line from file + read (14,'(A)',iostat=io_error) dummy + i = i + 1 + + ! if we're past the header AND the read worked (i.e. not end of file) + if (i > n_header .and. io_error==0) then + + k = count_seperators_in_line(dummy, sep) + if (k == ncols) then + dummy = adjustl(dummy) + read(dummy,*) array(i-n_header,:) + else + write(*,*) "ERROR in file ",trim(adjustl(file))," line ", i, "has wrong number of columns, skip" + ! this module is not MPI aware so abort is not available + stop + endif + endif + enddo + close (14) + + write(*,'("Done reading.")') + write(*,'(80("-"))') end subroutine read_array_from_ascii_file @@ -106,63 +125,96 @@ end subroutine read_array_from_ascii_file ! count the number of lines in an ascii file, skip n_header lines ! --> MPI wrapper in the MPI parser module !----------------------------------------------------------------------------- - subroutine count_lines_in_ascii_file(file, num_lines, n_header) - implicit none - character(len=*), intent(in) :: file - integer, intent(out) :: num_lines - integer, intent(in) :: n_header - integer :: io_error, i - character(len=maxcolumns) :: dummy - - ! check if the specified file exists - call check_file_exists( file ) + subroutine count_lines_in_ascii_file(file, num_lines, n_header, user_separator) + implicit none + character(len=*), intent(in) :: file + integer, intent(out) :: num_lines + integer, intent(in) :: n_header + character(len=1), optional, intent(in) :: user_separator + integer :: io_error, i + character(len=maxcolumns) :: dummy + character(len=1) :: sep + + ! check if the specified file exists + call check_file_exists( file ) + + if (present(user_separator)) then + sep = user_separator + else + sep = " " + endif - ! count the lines - io_error = 0 - i = 0 - open(unit=14,file=trim(adjustl(file)),action='read',status='old') - do while (io_error==0) - read (14,'(A)',iostat=io_error) dummy - if (io_error==0) i = i+1 - enddo - close (14) - num_lines = i - n_header + ! count the lines + io_error = 0 + i = 0 + open(unit=14,file=trim(adjustl(file)),action='read',status='old') + do while (io_error==0) + read (14,'(A)',iostat=io_error) dummy + if (io_error==0) i = i+1 + enddo + close (14) + num_lines = i - n_header end subroutine count_lines_in_ascii_file - subroutine count_cols_in_ascii_file(file, num_cols, n_header) - implicit none - character(len=*), intent(in) :: file - integer, intent(out) :: num_cols - integer, intent(in) :: n_header - integer :: io_error, i - character(len=maxcolumns) :: dummy + subroutine count_cols_in_ascii_file(file, num_cols, n_header, user_separator) + implicit none + character(len=*), intent(in) :: file + integer, intent(out) :: num_cols + integer, intent(in) :: n_header + character(len=1), optional, intent(in) :: user_separator + integer :: io_error, i + character(len=maxcolumns) :: dummy + character(len=1) :: sep - ! check if the specified file exists - call check_file_exists( file ) + ! check if the specified file exists + call check_file_exists( file ) - ! count the lines - io_error = 0 - i = 0 + if (present(user_separator)) then + sep = user_separator + else + sep = " " + endif - open(unit=14,file=trim(adjustl(file)),action='read',status='old') - do while (io_error==0) - read (14,'(A)',iostat=io_error) dummy - if (io_error==0 .and. i>n_header) exit - i=i+1 - enddo - close (14) + ! count the lines + io_error = 0 + i = 1 - num_cols = 1 - do i = 1, len_trim(adjustl(dummy)) - ! count elements in the line by counting the separating spaces - if ( dummy(i:i) == " " ) then - num_cols = num_cols + 1 - end if - enddo + ! just read first line....after the header + open(unit=14,file=trim(adjustl(file)),action='read',status='old') + do while (io_error==0) + read (14,'(A)',iostat=io_error) dummy + if (io_error==0 .and. i>n_header) exit + i = i+1 + enddo + close (14) + + num_cols = count_seperators_in_line(dummy, sep) end subroutine count_cols_in_ascii_file + + integer function count_seperators_in_line( line, sep ) + implicit none + CHARACTER(len=*), intent(in) :: line, sep + character(len=maxcolumns) :: dummy + integer :: num_cols, i + + num_cols = 1 + ! remove leading spaces + dummy = adjustl(line) + + do i = 1, len_trim(dummy) + ! count elements in the dummy by counting the separators + if ( dummy(i:i) == sep(1:1) ) then + num_cols = num_cols + 1 + end if + enddo + + count_seperators_in_line = num_cols + end function + + !------------------------------------------------------------------------------- ! clean a previously read ini file, deallocate its string array, and reset ! verbosity to .true. (as a matter of precaution) diff --git a/src/ini_files_parser_mpi.f90 b/src/ini_files_parser_mpi.f90 index 6ab69be..61e8e99 100644 --- a/src/ini_files_parser_mpi.f90 +++ b/src/ini_files_parser_mpi.f90 @@ -24,11 +24,12 @@ module module_ini_files_parser_mpi ! read an array from an ascii file (SERIAL version, to be executed only on root) ! note: array is assumed-shape and its size defines what we try to read !----------------------------------------------------------------------------- - subroutine read_array_from_ascii_file_mpi(file, array, n_header) + subroutine read_array_from_ascii_file_mpi(file, array, n_header, user_separator) implicit none character(len=*), intent(in) :: file integer, intent(in) :: n_header real(kind=pr), intent(inout) :: array (1:,1:) + character(len=1), optional, intent(in) :: user_separator integer :: nlines, ncols, mpicode, mpirank ! fetch my process id @@ -41,7 +42,12 @@ subroutine read_array_from_ascii_file_mpi(file, array, n_header) ncols = size(array,2) ! only root reads from file... - if (mpirank==0) call read_array_from_ascii_file(file, array, n_header) + if (present(user_separator)) then + if (mpirank==0) call read_array_from_ascii_file(file, array, n_header, user_separator) + else + if (mpirank==0) call read_array_from_ascii_file(file, array, n_header) + endif + ! ... then broadcast call MPI_BCAST(array,nlines*ncols,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpicode) end subroutine read_array_from_ascii_file_mpi @@ -50,35 +56,48 @@ end subroutine read_array_from_ascii_file_mpi !----------------------------------------------------------------------------- ! count the number of lines in an ascii file, skip n_header lines !----------------------------------------------------------------------------- - subroutine count_lines_in_ascii_file_mpi(file, num_lines, n_header) + subroutine count_lines_in_ascii_file_mpi(file, num_lines, n_header, user_separator) implicit none character(len=*), intent(in) :: file integer, intent(out) :: num_lines integer, intent(in) :: n_header + character(len=1), optional, intent(in) :: user_separator integer :: mpicode, mpirank ! fetch my process id call MPI_Comm_rank(MPI_COMM_WORLD, mpirank, mpicode) ! only root reads from file... - if (mpirank==0) call count_lines_in_ascii_file(file, num_lines, n_header) + if (present(user_separator)) then + if (mpirank==0) call count_lines_in_ascii_file(file, num_lines, n_header, user_separator) + else + if (mpirank==0) call count_lines_in_ascii_file(file, num_lines, n_header) + endif + ! ... then broadcast call MPI_BCAST(num_lines,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpicode) end subroutine count_lines_in_ascii_file_mpi - subroutine count_cols_in_ascii_file_mpi(file, num_cols, n_header) + subroutine count_cols_in_ascii_file_mpi(file, num_cols, n_header, user_separator) implicit none character(len=*), intent(in) :: file integer, intent(out) :: num_cols integer, intent(in) :: n_header + character(len=1), optional, intent(in) :: user_separator integer :: mpicode, mpirank ! fetch my process id call MPI_Comm_rank(MPI_COMM_WORLD, mpirank, mpicode) ! only root reads from file... - if (mpirank==0) call count_cols_in_ascii_file(file, num_cols, n_header) + if (present(user_separator)) then + if (mpirank==0) call count_cols_in_ascii_file(file, num_cols, n_header, user_separator) + else + if (mpirank==0) call count_cols_in_ascii_file(file, num_cols, n_header) + endif + + ! ... then broadcast call MPI_BCAST(num_cols,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpicode) end subroutine count_cols_in_ascii_file_mpi diff --git a/src/postprocessing/POD.f90 b/src/postprocessing/POD.f90 index edfb36d..408afee 100644 --- a/src/postprocessing/POD.f90 +++ b/src/postprocessing/POD.f90 @@ -1,8 +1,3 @@ -!------------------------------------------------------------------------------- -! ./flusi --postprocess --time-avg file_list.txt avgx_0000.h5 -! Reads in a list of files from a file, then loads one file after the other and -! computes the average field, which is then stored in the specified file. -!------------------------------------------------------------------------------- subroutine POD(help) use vars use p3dfft_wrapper @@ -12,24 +7,30 @@ subroutine POD(help) implicit none logical, intent(in) :: help character(len=strlen) :: fnamex_list, fname_this, fnamey_list, fnamez_list, dummy - character(len=strlen) :: fnamex, fnamey, fnamez - character(len=strlen) :: fname_list(1:3) + character(len=strlen) :: fnamex, fnamey, fnamez, timesteps_list + character(len=strlen) :: fname_list(1:3), fname_acoefs real(kind=pr), dimension(:,:,:), allocatable :: field_avg, field real(kind=pr), dimension(:,:), allocatable :: X_data, POD_modes, a_coefs - real(kind=pr), dimension(:), allocatable :: X_mean + real(kind=pr), dimension(:), allocatable :: X_mean, times DOUBLE PRECISION, dimension(:,:), allocatable :: C, V, D DOUBLE PRECISION, dimension(:), allocatable :: eigenvalues, work - integer :: ix, iy ,iz, io_error, i,j, N_modes, N_snapshots, info, it, dim, npoints - integer :: ncomponents + integer :: ix, iy ,iz, io_error, i,j, N_modes, N_snapshots, info, it, dim, npoints, k + integer :: ncomponents, nmodes_to_read, nmodes_file real(kind=pr) :: time, a, norm - LOGICAL :: vector + LOGICAL :: vector, reconstruct_all_time_steps, reconstruct_list, only_recon, total_not_fluctuations, only_modes if (help.and.root) then write(*,*) "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" - write(*,*) "./flusi -p --POD --components 3 --list file_list.txt [list_uy.txt] [list_uz.txt] --modes 10" + write(*,*) "./flusi -p --POD --components 3 --list file_list.txt [list_uy.txt] [list_uz.txt] " + write(*,*) "--modes 10 --reconstruct-all-time-steps --total-not-fluctuations" write(*,*) "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" write(*,*) " --list: a TXT file which simply contains the list of snapshots to read, one file per line" - write(*,*) " " + write(*,*) " --reconstruct-all-time-steps" + write(*,*) " --reconstruct-list timesteps.txt" + write(*,*) " --only-modes" + write(*,*) " --only-reconstruction a_coefs.txt (reconstruction from 34 modes, read from HDD)" + write(*,*) " --modes 10" + write(*,*) " --total-not-fluctuations" write(*,*) "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" write(*,*) "Parallel: Yes" return @@ -37,21 +38,50 @@ subroutine POD(help) ! defaults: ncomponents = 1 + reconstruct_all_time_steps = .false. + reconstruct_list = .false. + only_recon = .false. + total_not_fluctuations = .false. + only_modes = .false. ! fetch parameters from command line call do i = 1, COMMAND_ARGUMENT_COUNT() call get_command_argument(i,dummy) select case (dummy) + case ("--total-not-fluctuations") + total_not_fluctuations = .true. + if (root) write(*,*) "We do NOT remove the mean and work on TOTAL field" + + case ("--only-modes") + only_modes = .true. + if (root) write(*,*) "We compute only modes and do NOT reconstruct" + + case ("--only-reconstruction") + only_recon = .true. + + call get_command_argument(i+1, fname_acoefs) + if (root) write(*,*) "reconstruction from acoefs=", fname_acoefs + case ("--list") do j = 1, ncomponents call get_command_argument(i+j, fname_list(j)) enddo + case ("--reconstruct-all-time-steps") + reconstruct_all_time_steps = .true. + case ("--modes") call get_command_argument(i+1,dummy) read(dummy,*) N_modes if (root) write(*,*) "Will save N_modes=", N_modes + case ("--reconstruct-list") + call get_command_argument(i+1,timesteps_list) + reconstruct_list = .true. + reconstruct_all_time_steps = .false. + if (root) write(*,*) "Will reconstruct at time steps given by ", timesteps_list + call check_file_exists(timesteps_list) + case ("--components") call get_command_argument(i+1,dummy) read(dummy,*) ncomponents @@ -60,231 +90,372 @@ subroutine POD(help) end select enddo - !----------------------------------------------------------------------------- - ! check if input file exists, the file contains the list of h5 files to be avg - !----------------------------------------------------------------------------- - do j = 1, ncomponents - call check_file_exists ( fname_list(j) ) - if (root) write(*,*) "Reading list of files from "//fname_list(j) - enddo - + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! VARIANT A: read data, compute POD + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (.not. only_recon) then + !----------------------------------------------------------------------------- + ! check if input file exists, the file contains the list of h5 files to be avg + !----------------------------------------------------------------------------- + do j = 1, ncomponents + call check_file_exists ( fname_list(j) ) + if (root) write(*,*) "Reading list of files from "//fname_list(j) + enddo - !----------------------------------------------------------------------------- - ! read in the file, loop over lines - !----------------------------------------------------------------------------- - call count_lines_in_ascii_file_mpi(fname_list(1), N_snapshots, n_header=0) - if (root) write(*,*) "Reading N_snapshots=", N_snapshots + !----------------------------------------------------------------------------- + ! read in the file, loop over lines + !----------------------------------------------------------------------------- + call count_lines_in_ascii_file_mpi(fname_list(1), N_snapshots, n_header=0) - allocate( C(1:N_snapshots,1:N_snapshots) ) - allocate( D(1:N_snapshots,1:N_snapshots) ) - allocate( V(1:N_snapshots,1:N_snapshots) ) - allocate( eigenvalues(1:N_snapshots) ) - allocate( a_coefs(1:N_snapshots,1:N_snapshots) ) - allocate( work(1:5*N_snapshots) ) + if (root) write(*,*) "Reading N_snapshots=", N_snapshots - do j = 1, ncomponents - open( unit=10+j, file=fname_list(j), action='read', status='old' ) - enddo + allocate( C(1:N_snapshots,1:N_snapshots) ) + allocate( D(1:N_snapshots,1:N_snapshots) ) + allocate( V(1:N_snapshots,1:N_snapshots) ) + allocate( eigenvalues(1:N_snapshots) ) + allocate( a_coefs(1:N_snapshots,1:N_snapshots) ) + allocate( work(1:5*N_snapshots) ) + allocate( times(1:N_snapshots) ) - io_error = 0 - i = 1 - do while (i <= N_snapshots) do j = 1, ncomponents - read (10+j, '(A)', iostat=io_error) fnamex - call check_file_exists ( fnamex ) + open( unit=10+j, file=fname_list(j), action='read', status='old' ) + enddo - ! initialization is done after first read. - if ( i == 1 .and. j == 1 ) then - ! get file size etc - call fetch_attributes( fnamex, nx, ny, nz, xl, yl, zl, time, nu, origin ) - ! initialization parallel module (no FFTS) - call decomposition_initialize() + io_error = 0 + i = 1 + do while (i <= N_snapshots) + do j = 1, ncomponents + read (10+j, '(A)', iostat=io_error) fnamex + call check_file_exists ( fnamex ) - ! size of a flattened snapshot - npoints = (rb(1)-ra(1)+1) * (rb(2)-ra(2)+1) * (rb(3)-ra(3)+1) + ! initialization is done after first read. + if ( i == 1 .and. j == 1 ) then + ! get file size etc + call fetch_attributes( fnamex, nx, ny, nz, xl, yl, zl, time, nu, origin ) + ! initialization parallel module (no FFTS) + call decomposition_initialize() - if (decomposition /= "1D") call abort(28122018,"I think this module works only for 1D MPI decomp (scalar products)") + ! size of a flattened snapshot + npoints = (rb(1)-ra(1)+1) * (rb(2)-ra(2)+1) * (rb(3)-ra(3)+1) - ! memory for one field - allocate( field( ra(1):rb(1),ra(2):rb(2),ra(3):rb(3) ) ) + if (decomposition /= "1D") then + call abort(28122018,"I think this module works only for 1D MPI decomp (scalar products)") + endif - allocate( X_data( 1:npoints * ncomponents, 1:N_snapshots) ) - allocate( POD_modes( 1:npoints * ncomponents, 1:N_snapshots) ) - allocate( X_mean( 1:npoints * ncomponents ) ) - endif + ! memory for one field + allocate( field( ra(1):rb(1),ra(2):rb(2),ra(3):rb(3) ) ) - ! read the field from file - call read_single_file( fnamex, field ) + allocate( X_data( 1:npoints * ncomponents, 1:N_snapshots) ) + allocate( POD_modes( 1:npoints * ncomponents, 1:N_snapshots) ) + allocate( X_mean( 1:npoints * ncomponents ) ) + endif - ! add it to the data array - X_data( 1+(j-1)*npoints:(j)*npoints, i) = reshape(field, (/npoints/) ) - enddo - if (root) write(*,*) "filled snapshot slot", i - i = i+1 - enddo + ! read the field from file + call read_single_file( fnamex, field ) - do j = 1, ncomponents - close(10+j) - enddo + ! add it to the data array + X_data( 1+(j-1)*npoints:(j)*npoints, i) = reshape(field, (/npoints/) ) + enddo + if (root) write(*,*) "filled snapshot slot", i + i = i+1 + enddo - !--------------------------------------------------------------------------- - ! compute fluctuations (remove mean) - !--------------------------------------------------------------------------- + do j = 1, ncomponents + close(10+j) + enddo - ! compute mean of data - X_mean = 0.0_pr - do i = 1, N_snapshots - X_mean = X_mean + X_data(:,i) - enddo - X_mean = X_mean / dble(N_snapshots) + !--------------------------------------------------------------------------- + ! compute fluctuations (remove mean) + !--------------------------------------------------------------------------- + if ( .not. total_not_fluctuations ) then + ! compute mean of data + X_mean = 0.0_pr + do i = 1, N_snapshots + X_mean = X_mean + X_data(:,i) + enddo + X_mean = X_mean / dble(N_snapshots) - ! POD acts on fluctuations, so remove the mean. - do i = 1, N_snapshots - X_data(:,i) = X_data(:,i) - X_mean - enddo + ! POD acts on fluctuations, so remove the mean. + do i = 1, N_snapshots + X_data(:,i) = X_data(:,i) - X_mean + enddo + endif - ! divide by number of snapshots (eqn. 3.29 on p. 11 of Luchtenberg, Noak.) - X_data = X_data / dble(N_snapshots) + !--------------------------------------------------------------------------- + ! construction of covariance matrix from snapshot data + !--------------------------------------------------------------------------- + ! compute matrix C + C = matmul( transpose(X_data), X_data ) - !--------------------------------------------------------------------------- - ! construction of covariance matrix from snapsot data - !--------------------------------------------------------------------------- - ! compute matrix C - C = matmul( transpose(X_data), X_data ) - ! do i = 1, N_snapshots - ! do j = 1, N_snapshots - ! C(i,j) = sum( X_data(:,i)*X_data(:,j) ) - ! enddo - ! enddo - - ! normalization of C Matrix - C = C / dble( nx*ny*nz ) - call MPI_ALLREDUCE(MPI_IN_PLACE, C, N_snapshots**2, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) - - - ! if (root) then - ! open(14, file='C_matrix_fortran.txt', status='replace') - ! do i = 1, N_snapshots - ! write(14,'(400(es15.8,1x))') C(i,:) - ! enddo - ! close(14) - ! end if + ! divide by number of snapshots (eqn. 3.29 on p. 11 of Luchtenberg, Noak.) + C = C / dble(N_snapshots) - !--------------------------------------------------------------------------- - ! eigenvalues of covariance matrix - !--------------------------------------------------------------------------- - call DSYEV('V', 'U', N_snapshots, C, N_snapshots, eigenvalues, work, 5*N_snapshots, info) - ! as in matlab the eigenvalues are sorted in ascending order... - ! on output, C now contains the eigenvectors: - V = C + call MPI_ALLREDUCE(MPI_IN_PLACE, C, N_snapshots**2, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) - if (root) write(*,*) "info=", info - if (info /= 0) call abort(333,"The eigenvalue solver failed...") + ! normalization of C Matrix + C = C / dble( nx*ny*nz ) - if (root) then - write(*,*) "----v eigenvalues v-----" - do i = 1, N_snapshots - write(*,'(1(es12.4,1x))') eigenvalues(i) - enddo - write(*,*) "----^ eigenvalues ^-----" - endif + ! if (root) then + ! open(14, file='C_matrix_fortran.txt', status='replace') + ! do i = 1, N_snapshots + ! write(14,'(400(es15.8,1x))') C(i,:) + ! enddo + ! close(14) + ! end if - !--------------------------------------------------------------------------- - ! construct POD basis functions (modes) - !--------------------------------------------------------------------------- - ! eqn. 31 from AIAA review + !--------------------------------------------------------------------------- + ! eigenvalues of covariance matrix + !--------------------------------------------------------------------------- + call DSYEV('V', 'U', N_snapshots, C, N_snapshots, eigenvalues, work, 5*N_snapshots, info) + ! as in matlab the eigenvalues are sorted in ascending order... + ! on output, C now contains the eigenvectors: + V = C - POD_modes = matmul(X_data, V) + if (root) write(*,*) "info=", info + if (info /= 0) call abort(333,"The eigenvalue solver failed...") - ! do i = 1, size(POD_modes, 1) - ! do j = 1, N_snapshots - ! POD_modes(i,j) = sum( X_data(i,:)*V(:,j) ) - ! enddo - ! enddo - do i = 1, N_snapshots - ! note the division by sqrt(lambda) might cause numerical instabilities - ! if the eigenvalue is close to zero - if ( eigenvalues(i) > 1.0e-9_pr ) then - POD_modes(:,i) = POD_modes(:,i) / sqrt(eigenvalues(i)) + if (root) then + write(*,*) "----v eigenvalues v-----" + do i = 1, N_snapshots + write(*,'(1(es12.4,1x))') eigenvalues(i) + enddo + write(*,*) "----^ eigenvalues ^-----" endif - ! ! for some reason the modes are not normalized here, so we take care of that - ! ! I think the reason is the scalar product, which needs to be scaled by nx*ny*nz (and NOT npoints) - ! norm = sqrt( sum(POD_modes(:,i)**2) / dble(nx*ny*nz) ) - ! call MPI_ALLREDUCE(MPI_IN_PLACE, norm, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) - ! if(root) write(*,*) i, norm - ! if (norm > 1.0e-8_pr) then - ! POD_modes(:,i) = POD_modes(:,i) / norm - ! endif + !--------------------------------------------------------------------------- + ! construct POD basis functions (modes) + !--------------------------------------------------------------------------- + ! eqn. 31 from AIAA review - enddo + if (root) write(*,*) "Constructing POD modes (X*V)" + POD_modes = matmul(X_data, V) - !--------------------------------------------------------------------------- - ! save modes to disk - !--------------------------------------------------------------------------- - do i = 1, N_modes - do j = 1, ncomponents - field = reshape( POD_modes(1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1), & - (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) - ! create filename - write(fname_this, '("mode",i1,"_",i3.3,".h5")') j, i + if (root) write(*,*) "Constructing POD modes (sqrt(lambda))" + do i = 1, N_snapshots + ! note the division by sqrt(lambda) might cause numerical instabilities + ! if the eigenvalue is close to zero + if ( eigenvalues(i) > 1.0e-9_pr ) then + POD_modes(:,i) = POD_modes(:,i) / sqrt( dble(N_snapshots)*eigenvalues(i)) + ! POD_modes(:,i) = POD_modes(:,i) / sqrt(eigenvalues(i)) + endif - call save_field_hdf5 ( dble(i), fname_this, field ) + ! ! for some reason the modes are not normalized here, so we take care of that + ! ! I think the reason is the scalar product, which needs to be scaled by nx*ny*nz (and NOT npoints) + ! norm = sqrt( sum(POD_modes(:,i)**2) / dble(nx*ny*nz) ) + ! + ! call MPI_ALLREDUCE(MPI_IN_PLACE, norm, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) + ! + ! if (root) write(*,*) i, norm + ! + ! if (norm > 1.0e-8_pr) then + ! POD_modes(:,i) = POD_modes(:,i) / norm + ! endif enddo - enddo - !--------------------------------------------------------------------------- - ! temporal coefficients - !--------------------------------------------------------------------------- - a_coefs = 0.0_pr - do it = 1, N_snapshots + !--------------------------------------------------------------------------- + ! save modes to disk + !--------------------------------------------------------------------------- + if (root) write(*,*) "Saving POD modes to disk" + do i = 1, N_modes - ! scalar product (inner product) - a_coefs(it,i) = sum( X_data(:, it) * POD_modes(:, N_snapshots-i+1) ) + do j = 1, ncomponents + field = reshape( POD_modes(1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1), & + (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + + ! create filename + write(fname_this, '("mode",i1,"_",i3.3,".h5")') j, i + + call save_field_hdf5 ( dble(i), fname_this, field ) + enddo enddo - enddo - call MPI_ALLREDUCE(MPI_IN_PLACE, a_coefs, N_snapshots**2, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) - a_coefs = a_coefs / dble(nx*ny*nz)!*ncomponents) - ! a_coefs = a_coefs * dble(ncomponents) + !--------------------------------------------------------------------------- + ! temporal coefficients + !--------------------------------------------------------------------------- + if (root) write(*,*) "Computing temporal coefficients a" - if (root) then - write(*,*) 'writing temporal coefficients a_coefs.txt' - open(14, file='a_coefs.txt', status='replace') - do i = 1, N_snapshots - write(14,'(400(es15.8,1x))') a_coefs(i, 1:N_modes) + a_coefs = 0.0_pr + do it = 1, N_snapshots + do i = 1, N_modes + ! scalar product (inner product) + a_coefs(it,i) = sum( X_data(:, it) * POD_modes(:, N_snapshots-i+1) ) + enddo enddo - close(14) - end if + call MPI_ALLREDUCE(MPI_IN_PLACE, a_coefs, N_snapshots**2, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, info) + a_coefs = a_coefs / dble(nx*ny*nz) + + if (root) then + write(*,*) 'writing temporal coefficients a_coefs.txt' + open(14, file='a_coefs.txt', status='replace') + do i = 1, N_snapshots + write(14,'(400(es15.8,1x))') a_coefs(i, 1:N_modes) + enddo + close(14) + end if + + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! VARIANT B: read previously computed POD from files + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + else + + !**************** read-from-file-mode **************** + + call count_lines_in_ascii_file_mpi(fname_acoefs, N_snapshots, n_header=0) + call count_cols_in_ascii_file_mpi(fname_acoefs, nmodes_file, n_header=0) + + if (root) write(*,*) "a_coefs.txt:", N_snapshots, nmodes_file + + allocate( a_coefs(1:N_snapshots, 1:nmodes_file)) + call read_array_from_ascii_file_mpi(fname_acoefs, a_coefs, n_header=0) + + if (root) then + write(*,*) "~~~~v ", fname_acoefs + do k = 1, N_snapshots + write(*,'(400(es15.8,1x))') a_coefs(k, :) + enddo + write(*,*) "~~~~^ ", fname_acoefs + endif + + if (N_snapshots < N_modes) then + call abort(7771,"You want to reconstruct using more modes than the POD originally used.") + endif + + ! read modes, but only as many as we require for reconstruction + do i = 1, N_modes + do j = 1, ncomponents + ! create filename + write(fname_this, '("mode",i1,"_",i3.3,".h5")') j, i + + if (root) write(*,*) "reading mode: ", fname_this + + if (i==1 .and. j==1) then + ! get file size etc + call fetch_attributes( fname_this, nx, ny, nz, xl, yl, zl, time, nu, origin ) + ! initialization parallel module (no FFTS) + call decomposition_initialize() + + ! size of a flattened snapshot + npoints = (rb(1)-ra(1)+1) * (rb(2)-ra(2)+1) * (rb(3)-ra(3)+1) + + if (decomposition /= "1D") then + call abort(28122018,"I think this module works only for 1D MPI decomp (scalar products)") + endif + + ! memory for one field + allocate( field( ra(1):rb(1),ra(2):rb(2),ra(3):rb(3) ) ) + allocate( POD_modes( 1:npoints * ncomponents, 1:N_snapshots) ) + endif + + ! read the mode + call read_single_file( fname_this, field ) + + ! sort it in the array (note this step is not necessary but simplifies coding) + POD_modes( 1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1) = reshape(field, (/npoints/) ) + + enddo + enddo + + endif + + + ! at this point, we have the POD_modes and their temporal coefficients ready + + ! if not reconstructing, we're done now + if (only_modes) return + + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! RECONSTRUCTION + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + !--------------------------------------------------------------------------- ! reconstruction using N_modes !--------------------------------------------------------------------------- ! which time step to reconstruct at - it = 1 + if (reconstruct_all_time_steps .and. .not. reconstruct_list) then + do it = 1, N_snapshots + do j = 1, ncomponents + field = 0.0_pr - do j = 1, ncomponents - field = 0.0_pr + do i = 1, N_modes + field = field + a_coefs(it,i) * reshape( POD_modes(1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1), & + (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + enddo - do i = 1, N_modes - field = field + a_coefs(it,i) * reshape( POD_modes(1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1), & - (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + ! create filename + write(fname_this,'("reconstruction",i2.2,"-",i1,"_",i3.3,".h5")') N_modes, j, it + call save_field_hdf5 ( dble(it), fname_this, field ) + + if (.not. only_recon) then + ! for comparison, also save original data (which is of course only the fluctuating + ! part of the solution). can only be done when computing POD, not on reconstruction + field = reshape( X_data(1+(j-1)*npoints:(j)*npoints, it), & + (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + + write(fname_this,'("original",i2.2,"-",i1,"_",i3.3,".h5")') N_modes, j, it + call save_field_hdf5 ( dble(it), fname_this, field ) + endif + + enddo enddo - ! create filename - write(fname_this,'("reconstruction",i1,"_",i3.3,".h5")') j, N_modes - call save_field_hdf5 ( dble(it), fname_this, field ) + !--------------------------------------------------------------------------- + ! reconstruction of a list of time steps + !--------------------------------------------------------------------------- + elseif ( .not. reconstruct_all_time_steps .and. reconstruct_list ) then - ! for comparison, also save original data (which is of course only the fluctuating - ! part of the solution) - field = reshape( X_data(1+(j-1)*npoints:(j)*npoints, it), (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + io_error = 0 + open(17, file=timesteps_list, action='read', status='old' ) - write(fname_this,'("original",i1,"_",i3.3,".h5")') j, N_modes - call save_field_hdf5 ( dble(it), fname_this, field ) + do while (io_error == 0) + read(17, *, iostat=io_error) it - enddo + if (io_error == 0) then + + if (root) write(*,*) "reconstruction it=", it, "using", N_modes, ncomponents + + if (it<1 .or. it>size(a_coefs,1)) then + write(*,*) "it=", it, size(a_coefs,1) + call abort(8881,"You request reconstruction at an invalid time step.") + endif + + do j = 1, ncomponents + field = 0.0_pr + + do i = 1, N_modes + field = field + a_coefs(it,i) * reshape( POD_modes(1+(j-1)*npoints:(j)*npoints, N_snapshots-i+1), & + (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + enddo + + ! create filename + write(fname_this,'("reconstruction",i2.2,"-",i1,"_",i3.3,".h5")') N_modes, j, it + call save_field_hdf5 ( dble(it), fname_this, field ) + + if (.not. only_recon) then + ! for comparison, also save original data (which is of course only the fluctuating + ! part of the solution). can only be done when computing POD, not on reconstruction + field = reshape( X_data(1+(j-1)*npoints:(j)*npoints, it), & + (/(rb(1)-ra(1)+1), (rb(2)-ra(2)+1), (rb(3)-ra(3)+1)/) ) + + write(fname_this,'("original",i2.2,"-",i1,"_",i3.3,".h5")') N_modes, j, it + call save_field_hdf5 ( dble(it), fname_this, field ) + endif + enddo + endif + enddo + + close(17) + endif end subroutine POD diff --git a/src/postprocessing/postprocessing.f90 b/src/postprocessing/postprocessing.f90 index a863ae6..98f696d 100644 --- a/src/postprocessing/postprocessing.f90 +++ b/src/postprocessing/postprocessing.f90 @@ -37,6 +37,8 @@ subroutine postprocessing() ! check what to do !----------------- select case (postprocessing_mode) + case ("--force") + call post_force(help) case ("--POD") call POD(help) case ("--divergence") @@ -148,6 +150,7 @@ subroutine postprocessing() write(*,*) "--extract-subset" write(*,*) "--field-analysis" write(*,*) "--force-decomp" + write(*,*) "--force" write(*,*) "--flexible-wing-mask" write(*,*) "--gradient" write(*,*) "--hdf2bin"