From 22808785cf963ab21456e7c2765460839d017d22 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 14:04:27 -0400 Subject: [PATCH 01/22] added new cmfd solver that doesnt require Petsc --- src/DEPENDENCIES | 8 + src/cmfd_execute.F90 | 5 + src/cmfd_solver.F90 | 489 ++++++++++++++++++++++++++++++++ src/input_xml.F90 | 8 +- src/matrix_header.F90 | 69 ++++- src/utils/build_dependencies.py | 2 +- src/vector_header.F90 | 25 ++ 7 files changed, 594 insertions(+), 12 deletions(-) create mode 100644 src/cmfd_solver.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index eeb538fc16..dfbfe5920f 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -24,6 +24,7 @@ cmfd_data.o: tally_header.o cmfd_execute.o: cmfd_data.o cmfd_execute.o: cmfd_jfnk_solver.o cmfd_execute.o: cmfd_power_solver.o +cmfd_execute.o: cmfd_solver.o cmfd_execute.o: constants.o cmfd_execute.o: error.o cmfd_execute.o: global.o @@ -76,6 +77,13 @@ cmfd_slepc_solver.o: cmfd_prod_operator.o cmfd_slepc_solver.o: constants.o cmfd_slepc_solver.o: global.o +cmfd_solver.o: cmfd_loss_operator.o +cmfd_solver.o: cmfd_prod_operator.o +cmfd_solver.o: constants.o +cmfd_solver.o: global.o +cmfd_solver.o: matrix_header.o +cmfd_solver.o: vector_header.o + cross_section.o: ace_header.o cross_section.o: constants.o cross_section.o: error.o diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 17c8157f41..83257a870d 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -22,6 +22,7 @@ contains use cmfd_data, only: set_up_cmfd use cmfd_power_solver, only: cmfd_power_execute use cmfd_jfnk_solver, only: cmfd_jfnk_execute + use cmfd_solver, only: cmfd_solver_execute use error, only: warning, fatal_error ! CMFD single processor on master @@ -37,6 +38,7 @@ contains call process_cmfd_options() ! Call solver +#ifdef PETSC if (trim(cmfd_solver_type) == 'power') then call cmfd_power_execute() elseif (trim(cmfd_solver_type) == 'jfnk') then @@ -45,6 +47,9 @@ contains message = 'solver type became invalid after input processing' call fatal_error() end if +#else + call cmfd_solver_execute() +#endif ! Save k-effective cmfd % k_cmfd(current_batch) = cmfd % keff diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 new file mode 100644 index 0000000000..f8ee244017 --- /dev/null +++ b/src/cmfd_solver.F90 @@ -0,0 +1,489 @@ +module cmfd_solver + +! This module contains routines to execute the power iteration solver + + use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix + use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix + use matrix_header, only: Matrix + use vector_header, only: Vector + + implicit none + private + public :: cmfd_solver_execute + + logical :: iconv ! did the problem converged + real(8) :: k_n ! new k-eigenvalue + real(8) :: k_o ! old k-eigenvalue + real(8) :: ktol = 1.e-8_8 ! tolerance on keff + real(8) :: stol = 1.e-8_8 ! tolerance on source + real(8) :: norm_n ! current norm of source vector + real(8) :: norm_o ! old norm of source vector + real(8) :: kerr ! error in keff + real(8) :: serr ! error in source + logical :: adjoint_calc ! run an adjoint calculation + type(Matrix) :: loss ! cmfd loss matrix + type(Matrix) :: prod ! cmfd prod matrix + type(Vector) :: phi_n ! new flux vector + type(Vector) :: phi_o ! old flux vector + type(Vector) :: s_n ! new source vector + type(Vector) :: s_o ! old flux vector + type(Vector) :: serr_v ! error in source + +contains + +!=============================================================================== +! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD +!=============================================================================== + + subroutine cmfd_solver_execute(k_tol, s_tol, adjoint) + + use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve + + real(8), optional :: k_tol ! tolerance on keff + real(8), optional :: s_tol ! tolerance on source + logical, optional :: adjoint ! adjoint calc + + logical :: physical_adjoint = .false. + + ! Set tolerances if present + if (present(k_tol)) ktol = k_tol + if (present(s_tol)) stol = s_tol + + ! Check for adjoint execution + adjoint_calc = .false. + if (present(adjoint)) adjoint_calc = adjoint + + ! Check for physical adjoint + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') & + physical_adjoint = .true. + + ! Start timer for build + call time_cmfdbuild % start() + + ! Initialize matrices and vectors + call init_data(physical_adjoint) + + ! Check for mathematical adjoint calculation + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') & + call compute_adjoint() + + ! Stop timer for build + call time_cmfdbuild % stop() + + ! Begin power iteration + call time_cmfdsolve % start() + call execute_power_iter() + call time_cmfdsolve % stop() + + ! Extract results + call extract_results() + + ! Deallocate data + call finalize() + + end subroutine cmfd_solver_execute + +!=============================================================================== +! INIT_DATA allocates matrices and vectors for CMFD solution +!=============================================================================== + + subroutine init_data(adjoint) + + use constants, only: ONE, ZERO + use global, only: cmfd_write_matrices + + logical :: adjoint + + integer :: n ! problem size + real(8) :: guess ! initial guess + + ! Set up matrices + call init_loss_matrix(loss) + call init_prod_matrix(prod) + + ! Get problem size + n = loss % n + + ! Set up flux vectors + call phi_n % create(n) + call phi_o % create(n) + + ! Set up source vectors + call s_n % create(n) + call s_o % create(n) + call serr_v % create(n) + + ! Set initial guess + guess = ONE + phi_n % val = guess + phi_o % val = guess + k_n = guess + k_o = guess + + ! Fill in loss matrix + call build_loss_matrix(loss, adjoint=adjoint) + + ! Fill in production matrix + call build_prod_matrix(prod, adjoint=adjoint) + + ! Setup petsc for everything + call loss % assemble() + call prod % assemble() + call loss % write('loss.dat') + call prod % write('prod.dat') + + ! Set norms to 0 + norm_n = ZERO + norm_o = ZERO + + end subroutine init_data + +!=============================================================================== +! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem +!=============================================================================== + + subroutine compute_adjoint() + + use global, only: cmfd_write_matrices + + ! Transpose matrices + call loss % transpose() + call prod % transpose() + + ! Write out matrix in binary file (debugging) + if (cmfd_write_matrices) then + call loss % write_petsc_binary('adj_lossmat.bin') + call prod % write_petsc_binary('adj_prodmat.bin') + end if + + end subroutine compute_adjoint + +!=============================================================================== +! EXECUTE_POWER_ITER is the main power iteration routine +! for the cmfd calculation +!=============================================================================== + + subroutine execute_power_iter() + + integer :: i ! iteration counter + + ! Reset convergence flag + iconv = .false. + + ! Begin power iteration + do i = 1, 10000 + + ! Compute source vector + call prod % vector_multiply(phi_o, s_o) + + ! Normalize source vector + s_o % val = s_o % val / k_o + + ! Compute new flux vector + call cmfd_linsolver(s_o, phi_n) + + ! Compute new source vector + call prod % vector_multiply(phi_n, s_n) + + ! Compute new k-eigenvalue + k_n = sum(s_n % val) / sum(s_o % val) + + ! Renormalize the old source + s_o % val = s_o % val * k_o + + ! Check convergence + call convergence(i) + + ! Break loop if converged + if (iconv) exit + + ! Record old values + phi_o % val = phi_n % val + k_o = k_n + norm_o = norm_n + + end do + + end subroutine execute_power_iter + +!=============================================================================== +! CONVERGENCE checks the convergence of the CMFD problem +!=============================================================================== + + subroutine convergence(iter) + + use constants, only: ONE, TINY_BIT + use global, only: cmfd_power_monitor, master + use, intrinsic :: ISO_FORTRAN_ENV + + integer :: iter ! iteration number + + ! Reset convergence flag + iconv = .false. + + ! Calculate error in keff + kerr = abs(k_o - k_n)/k_n + + ! Calculate max error in source + where (s_n % val > TINY_BIT) + serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 + end where + serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) + + ! Check for convergence + if(kerr < ktol .and. serr < stol) iconv = .true. + + ! Save the L2 norm of the source + norm_n = serr + + ! Print out to user + if (cmfd_power_monitor .and. master) then + write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", & + &1PE12.5,T55, "src-error: ",1PE12.5)') iter, k_n, kerr, serr + end if + + end subroutine convergence + +!=============================================================================== +! CMFD_LINSOLVER solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver(b, x) + + use constants, only: ONE, ZERO + use global, only: cmfd + + type(Vector) :: b ! right hand side vector + type(Vector) :: x ! unknown vector + + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: matidx ! matrix index of row + integer :: d1idx ! index of row "1" diagonal + integer :: d2idx ! index of row "2" diagonal + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: icol ! iteration counter over columns + logical :: found ! did we find col + real(8) :: m11 ! block diagonal component 1,1 + real(8) :: m12 ! block diagonal component 1,2 + real(8) :: m21 ! block diagonal component 2,1 + real(8) :: m22 ! block diagonal component 2,2 + real(8) :: dm ! determinant of block diagonal + real(8) :: d11 ! inverse component 1,1 + real(8) :: d12 ! inverse component 1,2 + real(8) :: d21 ! inverse component 2,1 + real(8) :: d22 ! inverse component 2,2 + real(8) :: tmp1 ! temporary sum g1 + real(8) :: tmp2 ! temporary sum g2 + real(8) :: err ! error in convergence of solution + real(8) :: tol ! tolerance on final error + type(Vector) :: tmpx ! temporary solution vector + + ! Set tolerance + tol = 1.e-10_8 + + ! Dimensions + ng = 2 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = ng*nx*ny*nz + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ZLOOP: do k = 1, nz + + YLOOP: do j = 1, ny + + XLOOP: do i = 1, nx + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get starting row in matrix for this block + call indices_to_matrix(1, i, j, k, matidx, ng, nx, ny) + + ! Get the index of the diagonals for both rows + call loss % search_indices(matidx, matidx, d1idx, found) + + call loss % search_indices(matidx + 1, matidx + 1, d2idx, found) + + ! Get block diagonal + m11 = loss % val(d1idx) ! group 1 diagonal + m12 = loss % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) + m21 = loss % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) + m22 = loss % val(d2idx) ! group 2 diagonal + + ! Analytically invert the diagonal + dm = m11*m22 - m12*m21 + d11 = m22/dm + d12 = -m12/dm + d21 = -m21/dm + d22 = m11/dm + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + tmp2 = ZERO + do icol = loss % get_row(matidx), d1idx - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = loss % get_row(matidx + 1), d2idx - 2 + tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = d1idx + 2, loss % get_row(matidx + 1) - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = d2idx + 1, loss % get_row(matidx + 2) - 1 + tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + + ! Adjust with RHS vector + tmp1 = b % val(matidx) - tmp1 + tmp2 = b % val(matidx + 1) - tmp2 + + ! Solve for new x + x % val(matidx) = d11*tmp1 + d12*tmp2 + x % val(matidx + 1) = d21*tmp1 + d22*tmp2 + + end do XLOOP + + end do YLOOP + + end do ZLOOP + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + if (err < tol) exit + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver + +!=============================================================================== +! EXTRACT_RESULTS takes results and puts them in CMFD global data object +!=============================================================================== + + subroutine extract_results() + + use global, only: cmfd, cmfd_write_matrices, current_batch + + character(len=25) :: filename ! name of file to write data + integer :: n ! problem size + + ! Get problem size + n = loss % n + + ! Allocate in cmfd object if not already allocated + if (adjoint_calc) then + if (.not. allocated(cmfd%adj_phi)) allocate(cmfd%adj_phi(n)) + else + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + end if + + ! Save values + if (adjoint_calc) then + cmfd % adj_phi = phi_n % val + else + cmfd % phi = phi_n % val + end if + + ! Save eigenvalue + if(adjoint_calc) then + cmfd%adj_keff = k_n + else + cmfd%keff = k_n + end if + + ! Normalize phi to 1 + if (adjoint_calc) then + cmfd%adj_phi = cmfd%adj_phi/sqrt(sum(cmfd%adj_phi*cmfd%adj_phi)) + else + cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi)) + end if + + ! Save dominance ratio + cmfd % dom(current_batch) = norm_n/norm_o + + ! Write out results + if (cmfd_write_matrices) then + if (adjoint_calc) then + filename = 'adj_fluxvec.bin' + else + filename = 'fluxvec.bin' + end if + call phi_n % write_petsc_binary(filename) + end if + + end subroutine extract_results + +!=============================================================================== +! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix +!=============================================================================== + + subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny) + + use global, only: cmfd, cmfd_coremap + + integer :: matidx ! the index location in matrix + integer :: i ! current x index + integer :: j ! current y index + integer :: k ! current z index + integer :: g ! current group index + integer :: nx ! maximum number of x cells + integer :: ny ! maximum number of y cells + integer :: ng ! maximum number of groups + + ! Check if coremap is used + if (cmfd_coremap) then + + ! Get idx from core map + matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) + + else + + ! Compute index + matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1) + + end if + + end subroutine indices_to_matrix + +!=============================================================================== +! FINALIZE frees all memory associated with power iteration +!=============================================================================== + + subroutine finalize() + + ! Destroy all objects +#ifdef PETSC + call gmres % destroy() +#endif + call loss % destroy() + call prod % destroy() + call phi_n % destroy() + call phi_o % destroy() + call s_n % destroy() + call s_o % destroy() + call serr_v % destroy + + end subroutine finalize + +end module cmfd_solver diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f50a4557c5..13b82e9400 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -579,10 +579,10 @@ contains call lower_case(run_cmfd_) if (run_cmfd_ == 'true' .or. run_cmfd_ == '1') then cmfd_run = .true. -#ifndef PETSC - message = 'CMFD is not available, recompile OpenMC with PETSc' - call fatal_error() -#endif +!#ifndef PETSC +! message = 'CMFD is not available, recompile OpenMC with PETSc' +! call fatal_error() +!#endif end if end subroutine read_settings_xml diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 25f1df95a5..51308123f5 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -20,17 +20,19 @@ module matrix_header # endif logical :: petsc_active contains - procedure :: create => matrix_create - procedure :: destroy => matrix_destroy - procedure :: add_value => matrix_add_value - procedure :: new_row => matrix_new_row - procedure :: assemble => matrix_assemble - procedure :: get_row => matrix_get_row - procedure :: get_col => matrix_get_col + procedure :: create => matrix_create + procedure :: destroy => matrix_destroy + procedure :: add_value => matrix_add_value + procedure :: new_row => matrix_new_row + procedure :: assemble => matrix_assemble + procedure :: get_row => matrix_get_row + procedure :: get_col => matrix_get_col procedure :: setup_petsc => matrix_setup_petsc procedure :: write_petsc_binary => matrix_write_petsc_binary procedure :: transpose => matrix_transpose procedure :: vector_multiply => matrix_vector_multiply + procedure :: search_indices => matrix_search_indices + procedure :: write => matrix_write end type matrix integer :: petsc_err @@ -355,4 +357,57 @@ contains end subroutine matrix_vector_multiply +!=============================================================================== +! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row +!=============================================================================== + + subroutine matrix_search_indices(self, row, col, idx, found) + + class(Matrix) :: self + integer :: row + integer :: col + integer :: idx + logical :: found + + integer :: j + + found = .false. + + COLS: do j = self % get_row(row), self % get_row(row + 1) - 1 + + if (self % get_col(j) == col) then + idx = j + found = .true. + exit + end if + + end do COLS + + end subroutine matrix_search_indices + +!=============================================================================== +! MATRIX_WRITE writes a matrix to file +!=============================================================================== + + subroutine matrix_write(self, filename) + + character(*) :: filename + class(Matrix) :: self + + integer :: unit_ + integer :: i + integer :: j + + open(newunit=unit_, file=filename) + + do i = 1, self % n + do j = self % get_row(i), self % get_row(i + 1) - 1 + write(unit_,*) i, self % get_col(j), self % val(j) + end do + end do + + close(unit_) + + end subroutine matrix_write + end module matrix_header diff --git a/src/utils/build_dependencies.py b/src/utils/build_dependencies.py index 19ac978857..cefc59b937 100755 --- a/src/utils/build_dependencies.py +++ b/src/utils/build_dependencies.py @@ -12,7 +12,7 @@ for src in glob.iglob('*.F90'): open(src, 'r').read()) for name in d: if name in ['mpi', 'hdf5', 'h5lt', 'petscsys', 'petscmat', 'petscksp', - 'petscsnes', 'petscvec']: + 'petscsnes', 'petscvec', 'omp_lib']: continue if name.startswith('xml_data_'): name = name.replace('xml_data_', 'templates/') diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 47b865ebec..ed9ceda089 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -23,6 +23,7 @@ module vector_header procedure :: add_value => vector_add_value procedure :: setup_petsc => vector_setup_petsc procedure :: write_petsc_binary => vector_write_petsc_binary + procedure :: copy => vector_copy end type Vector integer :: petsc_err @@ -123,4 +124,28 @@ contains end subroutine vector_write_petsc_binary +!=============================================================================== +! VECTOR_COPY allocates a separate vector and copies +!=============================================================================== + + subroutine vector_copy(self, vectocopy) + + class(Vector), target :: self + type(Vector) :: vectocopy + + ! Preallocate vector + if (.not.allocated(self % data)) allocate(self % data(vectocopy % n)) + self % val => self % data(1:vectocopy % n) + + ! Set n + self % n = vectocopy % n + + ! Copy values + self % val = vectocopy % val + + ! Petsc is default not active + self % petsc_active = .false. + + end subroutine vector_copy + end module vector_header From c3446ef71450c8866fbc0e208141c903e7850e6e Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 15:00:22 -0400 Subject: [PATCH 02/22] added SOR --- src/cmfd_input.F90 | 3 +++ src/cmfd_solver.F90 | 20 ++++++++++++++++---- src/global.F90 | 3 +++ src/templates/cmfd_t.xml | 1 + 4 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 9a768f4e3b..5c7b7355a4 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -199,6 +199,9 @@ contains cmfd_display = '' end if + ! Read in spectral radius estimate + cmfd_spectral = spectral_ + ! Create tally objects call create_cmfd_tally() diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index f8ee244017..53b3f8aed4 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -251,7 +251,7 @@ contains subroutine cmfd_linsolver(b, x) use constants, only: ONE, ZERO - use global, only: cmfd + use global, only: cmfd, cmfd_spectral type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector @@ -282,12 +282,16 @@ contains real(8) :: d22 ! inverse component 2,2 real(8) :: tmp1 ! temporary sum g1 real(8) :: tmp2 ! temporary sum g2 + real(8) :: x1 ! new g1 value of x + real(8) :: x2 ! new g2 value of x real(8) :: err ! error in convergence of solution real(8) :: tol ! tolerance on final error + real(8) :: w ! overrelaxation parameter type(Vector) :: tmpx ! temporary solution vector - ! Set tolerance + ! Set tolerance and overrelaxation parameter tol = 1.e-10_8 + w = ONE ! Dimensions ng = 2 @@ -356,8 +360,12 @@ contains tmp2 = b % val(matidx + 1) - tmp2 ! Solve for new x - x % val(matidx) = d11*tmp1 + d12*tmp2 - x % val(matidx + 1) = d21*tmp1 + d22*tmp2 + x1 = d11*tmp1 + d12*tmp2 + x2 = d21*tmp1 + d22*tmp2 + + ! Perform overrelaxation + x % val(matidx) = (ONE - w)*x % val(matidx) + w*x1 + x % val(matidx + 1) = (ONE - w)*x % val(matidx + 1) + w*x2 end do XLOOP @@ -369,8 +377,12 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) +print *, igs, err, w if (err < tol) exit + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + end do GS call tmpx % destroy() diff --git a/src/global.F90 b/src/global.F90 index f12171b257..79cddbaa47 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -355,6 +355,9 @@ module global ! CMFD display info character(len=25) :: cmfd_display + ! Estimate of spectral radius of CMFD matrices + real(8) :: cmfd_spectral = ONE + ! Information about state points to be written integer :: n_state_points = 0 type(SetInt) :: statepoint_batch diff --git a/src/templates/cmfd_t.xml b/src/templates/cmfd_t.xml index 0e49e0e5ab..901b30181d 100644 --- a/src/templates/cmfd_t.xml +++ b/src/templates/cmfd_t.xml @@ -30,5 +30,6 @@ + From c479adc5678506e06f5d79ab3cfa0cddb8ffb96b Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 18:00:12 -0400 Subject: [PATCH 03/22] solver now works with coremap feature, loop is around rows now instead of i,j,k --- src/cmfd_solver.F90 | 144 +++++++++++++++++++++++++++----------------- 1 file changed, 88 insertions(+), 56 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 53b3f8aed4..698e1dd8c1 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -256,6 +256,7 @@ contains type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector + integer :: g ! group index integer :: i ! loop counter for x integer :: j ! loop counter for y integer :: k ! loop counter for z @@ -264,11 +265,11 @@ contains integer :: ny ! maximum dimension in y direction integer :: nz ! maximum dimension in z direction integer :: ng ! number of energy groups - integer :: matidx ! matrix index of row integer :: d1idx ! index of row "1" diagonal integer :: d2idx ! index of row "2" diagonal integer :: igs ! Gauss-Seidel iteration counter integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration integer :: icol ! iteration counter over columns logical :: found ! did we find col real(8) :: m11 ! block diagonal component 1,1 @@ -298,7 +299,7 @@ contains nx = cmfd % indices(1) ny = cmfd % indices(2) nz = cmfd % indices(3) - n = ng*nx*ny*nz + n = loss % n ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 @@ -309,75 +310,67 @@ contains ! Perform red/black gs iterations REDBLACK: do irb = 0,1 - ZLOOP: do k = 1, nz + ! Begin loop around matrix rows + ROWS: do irow = 1, n, 2 - YLOOP: do j = 1, ny + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) - XLOOP: do i = 1, nx + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle - ! Filter out black cells (even) - if (mod(i+j+k,2) == irb) cycle + ! Get the index of the diagonals for both rows + call loss % search_indices(irow, irow, d1idx, found) + call loss % search_indices(irow + 1, irow + 1, d2idx, found) - ! Get starting row in matrix for this block - call indices_to_matrix(1, i, j, k, matidx, ng, nx, ny) + ! Get block diagonal + m11 = loss % val(d1idx) ! group 1 diagonal + m12 = loss % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) + m21 = loss % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) + m22 = loss % val(d2idx) ! group 2 diagonal - ! Get the index of the diagonals for both rows - call loss % search_indices(matidx, matidx, d1idx, found) + ! Analytically invert the diagonal + dm = m11*m22 - m12*m21 + d11 = m22/dm + d12 = -m12/dm + d21 = -m21/dm + d22 = m11/dm - call loss % search_indices(matidx + 1, matidx + 1, d2idx, found) + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + tmp2 = ZERO + do icol = loss % get_row(irow), d1idx - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = loss % get_row(irow + 1), d2idx - 2 + tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = d1idx + 2, loss % get_row(irow + 1) - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = d2idx + 1, loss % get_row(irow + 2) - 1 + tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + end do - ! Get block diagonal - m11 = loss % val(d1idx) ! group 1 diagonal - m12 = loss % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) - m21 = loss % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) - m22 = loss % val(d2idx) ! group 2 diagonal + ! Adjust with RHS vector + tmp1 = b % val(irow) - tmp1 + tmp2 = b % val(irow + 1) - tmp2 - ! Analytically invert the diagonal - dm = m11*m22 - m12*m21 - d11 = m22/dm - d12 = -m12/dm - d21 = -m21/dm - d22 = m11/dm + ! Solve for new x + x1 = d11*tmp1 + d12*tmp2 + x2 = d21*tmp1 + d22*tmp2 - ! Perform temporary sums, first do left of diag block, then right of diag block - tmp1 = ZERO - tmp2 = ZERO - do icol = loss % get_row(matidx), d1idx - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) - end do - do icol = loss % get_row(matidx + 1), d2idx - 2 - tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) - end do - do icol = d1idx + 2, loss % get_row(matidx + 1) - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) - end do - do icol = d2idx + 1, loss % get_row(matidx + 2) - 1 - tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) - end do + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + x % val(irow + 1) = (ONE - w)*x % val(irow + 1) + w*x2 - ! Adjust with RHS vector - tmp1 = b % val(matidx) - tmp1 - tmp2 = b % val(matidx + 1) - tmp2 - - ! Solve for new x - x1 = d11*tmp1 + d12*tmp2 - x2 = d21*tmp1 + d22*tmp2 - - ! Perform overrelaxation - x % val(matidx) = (ONE - w)*x % val(matidx) + w*x1 - x % val(matidx + 1) = (ONE - w)*x % val(matidx + 1) + w*x2 - - end do XLOOP - - end do YLOOP - - end do ZLOOP + end do ROWS end do REDBLACK ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) -print *, igs, err, w +! print *, igs, err, w if (err < tol) exit ! Calculation new overrelaxation parameter @@ -478,6 +471,45 @@ print *, igs, err, w end subroutine indices_to_matrix +!=============================================================================== +! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies +!=============================================================================== + + subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + use global, only: cmfd, cmfd_coremap + + integer :: i ! iteration counter for x + integer :: j ! iteration counter for y + integer :: k ! iteration counter for z + integer :: g ! iteration counter for groups + integer :: irow ! iteration counter over row (0 reference) + integer :: nx ! maximum number of x cells + integer :: ny ! maximum number of y cells + integer :: nz ! maximum number of z cells + integer :: ng ! maximum number of groups + + ! Check for core map + if (cmfd_coremap) then + + ! Get indices from indexmap + g = mod(irow-1, ng) + 1 + i = cmfd % indexmap((irow-1)/ng+1,1) + j = cmfd % indexmap((irow-1)/ng+1,2) + k = cmfd % indexmap((irow-1)/ng+1,3) + + else + + ! Compute indices + g = mod(irow-1, ng) + 1 + i = mod(irow-1, ng*nx)/ng + 1 + j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1 + k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1 + + end if + + end subroutine matrix_to_indices + !=============================================================================== ! FINALIZE frees all memory associated with power iteration !=============================================================================== From da9b1051b3b8b49d4c0f9a3608c33d309a45463a Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 19:06:15 -0400 Subject: [PATCH 04/22] default spectral radius is 0 (G-S no SOR) --- src/global.F90 | 2 +- src/templates/cmfd_t.xml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 79cddbaa47..62b196319e 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -356,7 +356,7 @@ module global character(len=25) :: cmfd_display ! Estimate of spectral radius of CMFD matrices - real(8) :: cmfd_spectral = ONE + real(8) :: cmfd_spectral = ZERO ! Information about state points to be written integer :: n_state_points = 0 diff --git a/src/templates/cmfd_t.xml b/src/templates/cmfd_t.xml index 901b30181d..093edf112b 100644 --- a/src/templates/cmfd_t.xml +++ b/src/templates/cmfd_t.xml @@ -30,6 +30,6 @@ - + From e9b0642c1f8959cfbc78da103fbe46f1e7551e36 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 19:06:59 -0400 Subject: [PATCH 05/22] added 1g solver and procedure pointer --- src/cmfd_solver.F90 | 132 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 127 insertions(+), 5 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 698e1dd8c1..0769beef73 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -29,6 +29,16 @@ module cmfd_solver type(Vector) :: s_o ! old flux vector type(Vector) :: serr_v ! error in source + ! CMFD linear solver interface + procedure(linsolve), pointer :: cmfd_linsolver => null() + abstract interface + subroutine linsolve(b, x) + import :: Vector + type(Vector) :: b + type(Vector) :: x + end subroutine linsolve + end interface + contains !=============================================================================== @@ -90,7 +100,8 @@ contains subroutine init_data(adjoint) use constants, only: ONE, ZERO - use global, only: cmfd_write_matrices + use error, only: fatal_error + use global, only: cmfd, cmfd_write_matrices, message logical :: adjoint @@ -136,6 +147,18 @@ contains norm_n = ZERO norm_o = ZERO + ! Set up solver + select case(cmfd % indices(4)) + case(1) + cmfd_linsolver => cmfd_linsolver_1g + case(2) + cmfd_linsolver => cmfd_linsolver_2g + case default + message = 'Must use PETSc for more than 2 groups' + call fatal_error() + end select + + end subroutine init_data !=============================================================================== @@ -245,10 +268,109 @@ contains end subroutine convergence !=============================================================================== -! CMFD_LINSOLVER solves the CMFD linear system +! CMFD_LINSOLVER_1g solves the CMFD linear system !=============================================================================== - subroutine cmfd_linsolver(b, x) + subroutine cmfd_linsolver_1g(b, x) + + use constants, only: ONE, ZERO + use global, only: cmfd, cmfd_spectral + + type(Vector) :: b ! right hand side vector + type(Vector) :: x ! unknown vector + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: tol ! tolerance on final error + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set tolerance and overrelaxation parameter + tol = 1.e-10_8 + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = loss % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get the index of the diagonals for both rows + call loss % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = loss % get_row(irow), didx - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + do icol = didx + 1, loss % get_row(irow + 1) - 1 + tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/loss % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + end do ROWS + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) +!print *, igs, err, w + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_1g + +!=============================================================================== +! CMFD_LINSOLVER_2G solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_2g(b, x) use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral @@ -370,7 +492,7 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) -! print *, igs, err, w +!print *, igs, err, w if (err < tol) exit ! Calculation new overrelaxation parameter @@ -380,7 +502,7 @@ contains call tmpx % destroy() - end subroutine cmfd_linsolver + end subroutine cmfd_linsolver_2g !=============================================================================== ! EXTRACT_RESULTS takes results and puts them in CMFD global data object From 035534124763d6b5a2a5b21ad064649736ff3886 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 24 Sep 2013 19:18:42 -0400 Subject: [PATCH 06/22] can now set G-S tolerance, updated tests --- src/cmfd_input.F90 | 3 ++- src/cmfd_solver.F90 | 4 ++-- src/global.F90 | 3 ++- src/templates/cmfd_t.xml | 1 + tests/test_cmfd_feed/cmfd.xml | 1 + tests/test_cmfd_nofeed/cmfd.xml | 1 + 6 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 5c7b7355a4..386b06fc27 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -199,8 +199,9 @@ contains cmfd_display = '' end if - ! Read in spectral radius estimate + ! Read in spectral radius estimate and G-S tolerance cmfd_spectral = spectral_ + cmfd_gs_tol = gs_tol_ ! Create tally objects call create_cmfd_tally() diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 0769beef73..b1b95b9823 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -274,7 +274,7 @@ contains subroutine cmfd_linsolver_1g(b, x) use constants, only: ONE, ZERO - use global, only: cmfd, cmfd_spectral + use global, only: cmfd, cmfd_spectral, cmfd_gs_tol type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector @@ -302,7 +302,7 @@ contains type(Vector) :: tmpx ! temporary solution vector ! Set tolerance and overrelaxation parameter - tol = 1.e-10_8 + tol = cmfd_gs_tol w = ONE ! Dimensions diff --git a/src/global.F90 b/src/global.F90 index 62b196319e..83d99e1f51 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -355,8 +355,9 @@ module global ! CMFD display info character(len=25) :: cmfd_display - ! Estimate of spectral radius of CMFD matrices + ! Estimate of spectral radius of CMFD matrices and G-S tolerance real(8) :: cmfd_spectral = ZERO + real(8) :: cmfd_gs_tol = 1.e-10_8 ! Information about state points to be written integer :: n_state_points = 0 diff --git a/src/templates/cmfd_t.xml b/src/templates/cmfd_t.xml index 093edf112b..04326bfdc9 100644 --- a/src/templates/cmfd_t.xml +++ b/src/templates/cmfd_t.xml @@ -31,5 +31,6 @@ + diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml index 51def74088..891d264167 100644 --- a/tests/test_cmfd_feed/cmfd.xml +++ b/tests/test_cmfd_feed/cmfd.xml @@ -12,5 +12,6 @@ dominance power true + 1.e-15 diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml index c36a396043..710b3a7d16 100644 --- a/tests/test_cmfd_nofeed/cmfd.xml +++ b/tests/test_cmfd_nofeed/cmfd.xml @@ -12,5 +12,6 @@ dominance power false + 1.e-15 From 66b5ede5c451a91be058b08285e2d1aa0dc44192 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 25 Sep 2013 09:16:06 -0400 Subject: [PATCH 07/22] added Wielandt shifting --- src/cmfd_input.F90 | 1 + src/cmfd_solver.F90 | 143 ++++++++++++++++++++++++++++++++------------ src/global.F90 | 1 + 3 files changed, 107 insertions(+), 38 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 386b06fc27..444c479599 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -202,6 +202,7 @@ contains ! Read in spectral radius estimate and G-S tolerance cmfd_spectral = spectral_ cmfd_gs_tol = gs_tol_ + if (shift_ /= ZERO) cmfd_shift = shift_ ! Create tally objects call create_cmfd_tally() diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index b1b95b9823..2fc74df039 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -14,6 +14,9 @@ module cmfd_solver logical :: iconv ! did the problem converged real(8) :: k_n ! new k-eigenvalue real(8) :: k_o ! old k-eigenvalue + real(8) :: k_s ! shift of eigenvalue + real(8) :: k_ln ! new shifted eigenvalue + real(8) :: k_lo ! old shifted eigenvalue real(8) :: ktol = 1.e-8_8 ! tolerance on keff real(8) :: stol = 1.e-8_8 ! tolerance on source real(8) :: norm_n ! current norm of source vector @@ -32,10 +35,13 @@ module cmfd_solver ! CMFD linear solver interface procedure(linsolve), pointer :: cmfd_linsolver => null() abstract interface - subroutine linsolve(b, x) + subroutine linsolve(A, b, x, i) + import :: Matrix import :: Vector + type(Matrix) :: A type(Vector) :: b type(Vector) :: x + integer :: i end subroutine linsolve end interface @@ -101,12 +107,13 @@ contains use constants, only: ONE, ZERO use error, only: fatal_error - use global, only: cmfd, cmfd_write_matrices, message + use global, only: cmfd, cmfd_write_matrices, message, cmfd_shift, keff logical :: adjoint integer :: n ! problem size real(8) :: guess ! initial guess + real(8) :: dw ! eigenvalue shift ! Set up matrices call init_loss_matrix(loss) @@ -128,8 +135,12 @@ contains guess = ONE phi_n % val = guess phi_o % val = guess - k_n = guess - k_o = guess + k_n = keff + k_o = k_n + dw = cmfd_shift + k_s = k_o + dw + k_ln = ONE/(ONE/k_n - ONE/k_s) + k_lo = k_ln ! Fill in loss matrix call build_loss_matrix(loss, adjoint=adjoint) @@ -188,11 +199,20 @@ contains subroutine execute_power_iter() + use constants, only: ONE + integer :: i ! iteration counter + integer :: innerits ! # of inner iterations + integer :: totalits ! total number of inners + type(Matrix) :: Ms ! Reset convergence flag iconv = .false. + ! Perform shift + call wielandt_shift(Ms) + totalits = 0 + ! Begin power iteration do i = 1, 10000 @@ -200,22 +220,26 @@ contains call prod % vector_multiply(phi_o, s_o) ! Normalize source vector - s_o % val = s_o % val / k_o + s_o % val = s_o % val / k_lo ! Compute new flux vector - call cmfd_linsolver(s_o, phi_n) + call cmfd_linsolver(Ms, s_o, phi_n, innerits) ! Compute new source vector call prod % vector_multiply(phi_n, s_n) - ! Compute new k-eigenvalue - k_n = sum(s_n % val) / sum(s_o % val) + ! Compute new shifted eigenvalue + k_ln = sum(s_n % val) / sum(s_o % val) + + ! Compute new eigenvalue + k_n = ONE/(ONE/k_ln + ONE/k_s) ! Renormalize the old source - s_o % val = s_o % val * k_o + s_o % val = s_o % val * k_lo ! Check convergence - call convergence(i) + call convergence(i, innerits) + totalits = totalits + innerits ! Break loop if converged if (iconv) exit @@ -223,23 +247,58 @@ contains ! Record old values phi_o % val = phi_n % val k_o = k_n + k_lo = k_ln norm_o = norm_n - +read* end do +print *, 'TOTAL INNER:', totalits + ! Destroy shifted matrix + call Ms % destroy() end subroutine execute_power_iter +!=============================================================================== +! WIELANDT SHIFT +!=============================================================================== + + subroutine wielandt_shift(Ms) + + use constants, only: ONE + + type(Matrix) :: Ms + + integer :: irow ! row counter + integer :: icol ! col counter + integer :: jcol ! current col index in prod matrix + + ! copy loss matrix + call Ms % copy(loss) + + ! perform subtraction + jcol = 1 + ROWS: do irow = 1, Ms % n + COLS: do icol = Ms % get_row(irow), Ms % get_row(irow + 1) - 1 + if (Ms % get_col(icol) == prod % get_col(jcol)) then + Ms % val(icol) = Ms % val(icol) - ONE/k_s*prod % val(jcol) + jcol = jcol + 1 + end if + end do COLS + end do ROWS + + end subroutine wielandt_shift + !=============================================================================== ! CONVERGENCE checks the convergence of the CMFD problem !=============================================================================== - subroutine convergence(iter) + subroutine convergence(iter, innerits) use constants, only: ONE, TINY_BIT use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV - integer :: iter ! iteration number + integer :: iter ! outer iteration number + integer :: innerits ! inner iteration nubmer ! Reset convergence flag iconv = .false. @@ -262,7 +321,8 @@ contains ! Print out to user if (cmfd_power_monitor .and. master) then write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", & - &1PE12.5,T55, "src-error: ",1PE12.5)') iter, k_n, kerr, serr + &1PE12.5,T55, "src-error: ",1PE12.5,T80,I0)') iter, k_n, kerr, & + serr, innerits end if end subroutine convergence @@ -271,13 +331,15 @@ contains ! CMFD_LINSOLVER_1g solves the CMFD linear system !=============================================================================== - subroutine cmfd_linsolver_1g(b, x) + subroutine cmfd_linsolver_1g(A, b, x, its) use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral, cmfd_gs_tol + type(Matrix) :: A ! coefficient matrix type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector + integer :: its ! number of inner iterations integer :: g ! group index integer :: i ! loop counter for x @@ -310,7 +372,7 @@ contains nx = cmfd % indices(1) ny = cmfd % indices(2) nz = cmfd % indices(3) - n = loss % n + n = A % n ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 @@ -331,19 +393,19 @@ contains if (mod(i+j+k,2) == irb) cycle ! Get the index of the diagonals for both rows - call loss % search_indices(irow, irow, didx, found) + call A % search_indices(irow, irow, didx, found) ! Perform temporary sums, first do left of diag block, then right of diag block tmp1 = ZERO - do icol = loss % get_row(irow), didx - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) end do - do icol = didx + 1, loss % get_row(irow + 1) - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) end do ! Solve for new x - x1 = (b % val(irow) - tmp1)/loss % val(didx) + x1 = (b % val(irow) - tmp1)/A % val(didx) ! Perform overrelaxation x % val(irow) = (ONE - w)*x % val(irow) + w*x1 @@ -354,6 +416,7 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs !print *, igs, err, w if (err < tol) exit @@ -370,13 +433,15 @@ contains ! CMFD_LINSOLVER_2G solves the CMFD linear system !=============================================================================== - subroutine cmfd_linsolver_2g(b, x) + subroutine cmfd_linsolver_2g(A, b, x, its) use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral + type(Matrix) :: A ! coefficient matrix type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector + integer :: its ! number of inner iterations integer :: g ! group index integer :: i ! loop counter for x @@ -421,7 +486,7 @@ contains nx = cmfd % indices(1) ny = cmfd % indices(2) nz = cmfd % indices(3) - n = loss % n + n = A % n ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 @@ -442,14 +507,14 @@ contains if (mod(i+j+k,2) == irb) cycle ! Get the index of the diagonals for both rows - call loss % search_indices(irow, irow, d1idx, found) - call loss % search_indices(irow + 1, irow + 1, d2idx, found) + call A % search_indices(irow, irow, d1idx, found) + call A % search_indices(irow + 1, irow + 1, d2idx, found) ! Get block diagonal - m11 = loss % val(d1idx) ! group 1 diagonal - m12 = loss % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) - m21 = loss % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) - m22 = loss % val(d2idx) ! group 2 diagonal + m11 = A % val(d1idx) ! group 1 diagonal + m12 = A % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) + m21 = A % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) + m22 = A % val(d2idx) ! group 2 diagonal ! Analytically invert the diagonal dm = m11*m22 - m12*m21 @@ -461,17 +526,17 @@ contains ! Perform temporary sums, first do left of diag block, then right of diag block tmp1 = ZERO tmp2 = ZERO - do icol = loss % get_row(irow), d1idx - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = A % get_row(irow), d1idx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) end do - do icol = loss % get_row(irow + 1), d2idx - 2 - tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = A % get_row(irow + 1), d2idx - 2 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) end do - do icol = d1idx + 2, loss % get_row(irow + 1) - 1 - tmp1 = tmp1 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = d1idx + 2, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) end do - do icol = d2idx + 1, loss % get_row(irow + 2) - 1 - tmp2 = tmp2 + loss % val(icol)*x % val(loss % get_col(icol)) + do icol = d2idx + 1, A % get_row(irow + 2) - 1 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) end do ! Adjust with RHS vector @@ -492,7 +557,9 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs !print *, igs, err, w +!read * if (err < tol) exit ! Calculation new overrelaxation parameter diff --git a/src/global.F90 b/src/global.F90 index 83d99e1f51..5f7011b900 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -358,6 +358,7 @@ module global ! Estimate of spectral radius of CMFD matrices and G-S tolerance real(8) :: cmfd_spectral = ZERO real(8) :: cmfd_gs_tol = 1.e-10_8 + real(8) :: cmfd_shift = INFINITY ! Information about state points to be written integer :: n_state_points = 0 From 9123098acc4ef0fbc39d542cf2ca08a209e75091 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 25 Sep 2013 09:16:41 -0400 Subject: [PATCH 08/22] added matrix copy and option to read in shift --- src/matrix_header.F90 | 31 +++++++++++++++++++++++++++++++ src/templates/cmfd_t.xml | 1 + 2 files changed, 32 insertions(+) diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 51308123f5..40f515c28d 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -33,6 +33,7 @@ module matrix_header procedure :: vector_multiply => matrix_vector_multiply procedure :: search_indices => matrix_search_indices procedure :: write => matrix_write + procedure :: copy => matrix_copy end type matrix integer :: petsc_err @@ -410,4 +411,34 @@ contains end subroutine matrix_write +!=============================================================================== +! MATRIX_COPY copies a matrix +!=============================================================================== + + subroutine matrix_copy(self, mattocopy) + + class(Matrix) :: self + type(Matrix) :: mattocopy + + ! Set n and nnz + self % n_count = mattocopy % n_count + self % nz_count = mattocopy % nz_count + self % n = mattocopy % n + self % nnz = mattocopy % nnz + + ! Allocate vectors + if (.not.allocated(self % row)) allocate(self % row(self % n + 1)) + if (.not.allocated(self % col)) allocate(self % col(self % nnz)) + if (.not.allocated(self % val)) allocate(self % val(self % nnz)) + + ! Set PETSc active to false + self % petsc_active = .false. + + ! Copy over data + self % row = mattocopy % row + self % col = mattocopy % col + self % val = mattocopy % val + + end subroutine matrix_copy + end module matrix_header diff --git a/src/templates/cmfd_t.xml b/src/templates/cmfd_t.xml index 04326bfdc9..7003973679 100644 --- a/src/templates/cmfd_t.xml +++ b/src/templates/cmfd_t.xml @@ -32,5 +32,6 @@ + From f8fe3adf3e2f04aa666d98f62c251c24d9a361c6 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 25 Sep 2013 10:14:14 -0400 Subject: [PATCH 09/22] added relative tolerancing and now can be set on input --- src/cmfd_input.F90 | 7 +++-- src/cmfd_solver.F90 | 57 +++++++++++++++++++++++----------------- src/global.F90 | 7 +++-- src/templates/cmfd_t.xml | 5 +++- 4 files changed, 47 insertions(+), 29 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 444c479599..2b6154b026 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -199,10 +199,13 @@ contains cmfd_display = '' end if - ! Read in spectral radius estimate and G-S tolerance + ! Read in spectral radius estimate and tolerances cmfd_spectral = spectral_ - cmfd_gs_tol = gs_tol_ if (shift_ /= ZERO) cmfd_shift = shift_ + cmfd_ktol = ktol_ + cmfd_stol = stol_ + cmfd_atoli = atoli_ + cmfd_rtoli = rtoli_ ! Create tally objects call create_cmfd_tally() diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 2fc74df039..972f3bacee 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -17,12 +17,12 @@ module cmfd_solver real(8) :: k_s ! shift of eigenvalue real(8) :: k_ln ! new shifted eigenvalue real(8) :: k_lo ! old shifted eigenvalue - real(8) :: ktol = 1.e-8_8 ! tolerance on keff - real(8) :: stol = 1.e-8_8 ! tolerance on source real(8) :: norm_n ! current norm of source vector real(8) :: norm_o ! old norm of source vector real(8) :: kerr ! error in keff real(8) :: serr ! error in source + real(8) :: ktol ! tolerance on keff + real(8) :: stol ! tolerance on source logical :: adjoint_calc ! run an adjoint calculation type(Matrix) :: loss ! cmfd loss matrix type(Matrix) :: prod ! cmfd prod matrix @@ -35,12 +35,13 @@ module cmfd_solver ! CMFD linear solver interface procedure(linsolve), pointer :: cmfd_linsolver => null() abstract interface - subroutine linsolve(A, b, x, i) + subroutine linsolve(A, b, x, tol, i) import :: Matrix import :: Vector type(Matrix) :: A type(Vector) :: b type(Vector) :: x + real(8) :: tol integer :: i end subroutine linsolve end interface @@ -51,20 +52,14 @@ contains ! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD !=============================================================================== - subroutine cmfd_solver_execute(k_tol, s_tol, adjoint) + subroutine cmfd_solver_execute(adjoint) use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve - real(8), optional :: k_tol ! tolerance on keff - real(8), optional :: s_tol ! tolerance on source logical, optional :: adjoint ! adjoint calc logical :: physical_adjoint = .false. - ! Set tolerances if present - if (present(k_tol)) ktol = k_tol - if (present(s_tol)) stol = s_tol - ! Check for adjoint execution adjoint_calc = .false. if (present(adjoint)) adjoint_calc = adjoint @@ -107,7 +102,8 @@ contains use constants, only: ONE, ZERO use error, only: fatal_error - use global, only: cmfd, cmfd_write_matrices, message, cmfd_shift, keff + use global, only: cmfd, cmfd_write_matrices, message, cmfd_shift, keff, & + cmfd_ktol, cmfd_stol logical :: adjoint @@ -167,8 +163,11 @@ contains case default message = 'Must use PETSc for more than 2 groups' call fatal_error() - end select - + end select + + ! Set tolerances + ktol = cmfd_ktol + stol = cmfd_stol end subroutine init_data @@ -200,15 +199,24 @@ contains subroutine execute_power_iter() use constants, only: ONE + use global, only: cmfd_atoli, cmfd_rtoli integer :: i ! iteration counter integer :: innerits ! # of inner iterations integer :: totalits ! total number of inners + real(8) :: atoli ! absolute minimum tolerance + real(8) :: rtoli ! relative tolerance based on source conv + real(8) :: toli ! the current tolerance of inners type(Matrix) :: Ms ! Reset convergence flag iconv = .false. + ! Set up tolerances + atoli = cmfd_atoli + rtoli = cmfd_rtoli + toli = rtoli*100._8 + ! Perform shift call wielandt_shift(Ms) totalits = 0 @@ -223,7 +231,7 @@ contains s_o % val = s_o % val / k_lo ! Compute new flux vector - call cmfd_linsolver(Ms, s_o, phi_n, innerits) + call cmfd_linsolver(Ms, s_o, phi_n, toli, innerits) ! Compute new source vector call prod % vector_multiply(phi_n, s_n) @@ -249,9 +257,12 @@ contains k_o = k_n k_lo = k_ln norm_o = norm_n -read* + + ! Get new tolerance for inners + toli = max(atoli, rtoli*serr) + end do -print *, 'TOTAL INNER:', totalits + ! Destroy shifted matrix call Ms % destroy() @@ -331,14 +342,15 @@ print *, 'TOTAL INNER:', totalits ! CMFD_LINSOLVER_1g solves the CMFD linear system !=============================================================================== - subroutine cmfd_linsolver_1g(A, b, x, its) + subroutine cmfd_linsolver_1g(A, b, x, tol, its) use constants, only: ONE, ZERO - use global, only: cmfd, cmfd_spectral, cmfd_gs_tol + use global, only: cmfd, cmfd_spectral type(Matrix) :: A ! coefficient matrix type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector + real(8) :: tol ! tolerance on final error integer :: its ! number of inner iterations integer :: g ! group index @@ -359,12 +371,10 @@ print *, 'TOTAL INNER:', totalits real(8) :: tmp1 ! temporary sum g1 real(8) :: x1 ! new g1 value of x real(8) :: err ! error in convergence of solution - real(8) :: tol ! tolerance on final error real(8) :: w ! overrelaxation parameter type(Vector) :: tmpx ! temporary solution vector - ! Set tolerance and overrelaxation parameter - tol = cmfd_gs_tol + ! Set overrelaxation parameter w = ONE ! Dimensions @@ -433,7 +443,7 @@ print *, 'TOTAL INNER:', totalits ! CMFD_LINSOLVER_2G solves the CMFD linear system !=============================================================================== - subroutine cmfd_linsolver_2g(A, b, x, its) + subroutine cmfd_linsolver_2g(A, b, x, tol, its) use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral @@ -441,6 +451,7 @@ print *, 'TOTAL INNER:', totalits type(Matrix) :: A ! coefficient matrix type(Vector) :: b ! right hand side vector type(Vector) :: x ! unknown vector + real(8) :: tol ! tolerance on final error integer :: its ! number of inner iterations integer :: g ! group index @@ -473,12 +484,10 @@ print *, 'TOTAL INNER:', totalits real(8) :: x1 ! new g1 value of x real(8) :: x2 ! new g2 value of x real(8) :: err ! error in convergence of solution - real(8) :: tol ! tolerance on final error real(8) :: w ! overrelaxation parameter type(Vector) :: tmpx ! temporary solution vector ! Set tolerance and overrelaxation parameter - tol = 1.e-10_8 w = ONE ! Dimensions diff --git a/src/global.F90 b/src/global.F90 index 5f7011b900..5f0f4f2e6b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -355,10 +355,13 @@ module global ! CMFD display info character(len=25) :: cmfd_display - ! Estimate of spectral radius of CMFD matrices and G-S tolerance + ! Estimate of spectral radius of CMFD matrices and tolerances real(8) :: cmfd_spectral = ZERO - real(8) :: cmfd_gs_tol = 1.e-10_8 real(8) :: cmfd_shift = INFINITY + real(8) :: cmfd_ktol = 1.e-8_8 + real(8) :: cmfd_stol = 1.e-8_8 + real(8) :: cmfd_atoli = 1.e-10_8 + real(8) :: cmfd_rtoli = 1.e-5_8 ! Information about state points to be written integer :: n_state_points = 0 diff --git a/src/templates/cmfd_t.xml b/src/templates/cmfd_t.xml index 7003973679..f372e25b7c 100644 --- a/src/templates/cmfd_t.xml +++ b/src/templates/cmfd_t.xml @@ -31,7 +31,10 @@ - + + + + From ccfc73e42232b2100112c952224142af6c0c208d Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 25 Sep 2013 10:29:34 -0400 Subject: [PATCH 10/22] removed copy of loss matrix for shift because we can overwrite loss matrix --- src/cmfd_solver.F90 | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 972f3bacee..edcbaf8db2 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -207,7 +207,6 @@ contains real(8) :: atoli ! absolute minimum tolerance real(8) :: rtoli ! relative tolerance based on source conv real(8) :: toli ! the current tolerance of inners - type(Matrix) :: Ms ! Reset convergence flag iconv = .false. @@ -218,7 +217,7 @@ contains toli = rtoli*100._8 ! Perform shift - call wielandt_shift(Ms) + call wielandt_shift() totalits = 0 ! Begin power iteration @@ -231,7 +230,7 @@ contains s_o % val = s_o % val / k_lo ! Compute new flux vector - call cmfd_linsolver(Ms, s_o, phi_n, toli, innerits) + call cmfd_linsolver(loss, s_o, phi_n, toli, innerits) ! Compute new source vector call prod % vector_multiply(phi_n, s_n) @@ -263,34 +262,26 @@ contains end do - ! Destroy shifted matrix - call Ms % destroy() - end subroutine execute_power_iter !=============================================================================== ! WIELANDT SHIFT !=============================================================================== - subroutine wielandt_shift(Ms) + subroutine wielandt_shift() use constants, only: ONE - type(Matrix) :: Ms - integer :: irow ! row counter integer :: icol ! col counter integer :: jcol ! current col index in prod matrix - ! copy loss matrix - call Ms % copy(loss) - ! perform subtraction jcol = 1 - ROWS: do irow = 1, Ms % n - COLS: do icol = Ms % get_row(irow), Ms % get_row(irow + 1) - 1 - if (Ms % get_col(icol) == prod % get_col(jcol)) then - Ms % val(icol) = Ms % val(icol) - ONE/k_s*prod % val(jcol) + ROWS: do irow = 1, loss % n + COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1 + if (loss % get_col(icol) == prod % get_col(jcol)) then + loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol) jcol = jcol + 1 end if end do COLS From 61f94289e6cd42e6349e87c7dc717cd4a610fead Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 25 Sep 2013 12:56:50 -0400 Subject: [PATCH 11/22] updated cmfd test cases for tighter tolerances --- tests/test_cmfd_feed/cmfd.xml | 3 ++- tests/test_cmfd_nofeed/cmfd.xml | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml index 891d264167..390e558639 100644 --- a/tests/test_cmfd_feed/cmfd.xml +++ b/tests/test_cmfd_feed/cmfd.xml @@ -12,6 +12,7 @@ dominance power true - 1.e-15 + 1.e-15 + 1.e-20 diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml index 710b3a7d16..4a4bf30108 100644 --- a/tests/test_cmfd_nofeed/cmfd.xml +++ b/tests/test_cmfd_nofeed/cmfd.xml @@ -12,6 +12,7 @@ dominance power false - 1.e-15 + 1.e-15 + 1.e-20 From d09a98ca131abd74c05fb6aeca4dacfb85e192d3 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 27 Sep 2013 16:55:39 -0400 Subject: [PATCH 12/22] cmfd shift cant be infinity due to FPE, making it very large and fixed matrix subtraction in solver --- src/cmfd_solver.F90 | 3 ++- src/global.F90 | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index edcbaf8db2..614362e40d 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -280,7 +280,8 @@ contains jcol = 1 ROWS: do irow = 1, loss % n COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1 - if (loss % get_col(icol) == prod % get_col(jcol)) then + if (loss % get_col(icol) == prod % get_col(jcol) .and. & + jcol < prod % get_row(irow + 1)) then loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol) jcol = jcol + 1 end if diff --git a/src/global.F90 b/src/global.F90 index 5f0f4f2e6b..46f7ba562b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -357,7 +357,7 @@ module global ! Estimate of spectral radius of CMFD matrices and tolerances real(8) :: cmfd_spectral = ZERO - real(8) :: cmfd_shift = INFINITY + real(8) :: cmfd_shift = 1.e6 real(8) :: cmfd_ktol = 1.e-8_8 real(8) :: cmfd_stol = 1.e-8_8 real(8) :: cmfd_atoli = 1.e-10_8 From 02dc33dd7d6bae0b87fceaa6f3e46e39a5e5ef13 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 09:52:37 -0400 Subject: [PATCH 13/22] updated preprocessing directives --- src/cmfd_solver.F90 | 46 ++++++++++++------------------------------- src/matrix_header.F90 | 8 +++++--- src/vector_header.F90 | 4 +++- 3 files changed, 21 insertions(+), 37 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 614362e40d..cc27b98574 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -102,8 +102,11 @@ contains use constants, only: ONE, ZERO use error, only: fatal_error - use global, only: cmfd, cmfd_write_matrices, message, cmfd_shift, keff, & + use global, only: cmfd, message, cmfd_shift, keff, & cmfd_ktol, cmfd_stol +#ifdef PETSC + use global, only: cmfd_write_matrices +#endif logical :: adjoint @@ -177,6 +180,9 @@ contains subroutine compute_adjoint() + use error, only: fatal_error + use global, only: message +#ifdef PETSC use global, only: cmfd_write_matrices ! Transpose matrices @@ -188,6 +194,10 @@ contains call loss % write_petsc_binary('adj_lossmat.bin') call prod % write_petsc_binary('adj_prodmat.bin') end if +#else + message = 'Adjoint calculations only allowed with PETSc' + call fatal_error() +#endif end subroutine compute_adjoint @@ -624,43 +634,13 @@ contains else filename = 'fluxvec.bin' end if +#ifdef PETSC call phi_n % write_petsc_binary(filename) +#endif end if end subroutine extract_results -!=============================================================================== -! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix -!=============================================================================== - - subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny) - - use global, only: cmfd, cmfd_coremap - - integer :: matidx ! the index location in matrix - integer :: i ! current x index - integer :: j ! current y index - integer :: k ! current z index - integer :: g ! current group index - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: ng ! maximum number of groups - - ! Check if coremap is used - if (cmfd_coremap) then - - ! Get idx from core map - matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) - - else - - ! Compute index - matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1) - - end if - - end subroutine indices_to_matrix - !=============================================================================== ! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies !=============================================================================== diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index eed6a91fb9..3f5aa204f7 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -27,13 +27,15 @@ module matrix_header procedure :: assemble => matrix_assemble procedure :: get_row => matrix_get_row procedure :: get_col => matrix_get_col - procedure :: setup_petsc => matrix_setup_petsc - procedure :: write_petsc_binary => matrix_write_petsc_binary - procedure :: transpose => matrix_transpose procedure :: vector_multiply => matrix_vector_multiply procedure :: search_indices => matrix_search_indices procedure :: write => matrix_write procedure :: copy => matrix_copy +# ifdef PETSC + procedure :: setup_petsc => matrix_setup_petsc + procedure :: write_petsc_binary => matrix_write_petsc_binary + procedure :: transpose => matrix_transpose +# endif end type matrix #ifdef PETSC diff --git a/src/vector_header.F90 b/src/vector_header.F90 index aab9b2bde0..4d8cb2cf5c 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -21,9 +21,11 @@ module vector_header procedure :: create => vector_create procedure :: destroy => vector_destroy procedure :: add_value => vector_add_value + procedure :: copy => vector_copy +# ifdef PETSC procedure :: setup_petsc => vector_setup_petsc procedure :: write_petsc_binary => vector_write_petsc_binary - procedure :: copy => vector_copy +# endif end type Vector integer :: petsc_err ! petsc error code From cf8c51d5af0f845cad2a0b38d94cd9c4ff2e95ad Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 10:02:44 -0400 Subject: [PATCH 14/22] removed some comments and fixed compilation error --- src/cmfd_solver.F90 | 3 --- src/input_xml.F90 | 6 ------ 2 files changed, 9 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index cc27b98574..49efe988dd 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -687,9 +687,6 @@ contains subroutine finalize() ! Destroy all objects -#ifdef PETSC - call gmres % destroy() -#endif call loss % destroy() call prod % destroy() call phi_n % destroy() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e615f9c16e..039e0419c3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -777,12 +777,6 @@ contains call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then cmfd_run = .true. -!#ifndef PETSC -! if (master) then -! message = 'CMFD is not available, compile OpenMC with PETSc' -! call fatal_error() -! end if -!#endif end if end if From 7398a1e8004b105b47f6d13f463e284f62d73b93 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 10:07:59 -0400 Subject: [PATCH 15/22] fixed misspelling for cmfd_matrices --- src/cmfd_input.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 09f02e61b3..61d6f6eaba 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -190,7 +190,7 @@ contains ! Output logicals if (check_for_node(doc, "write_matrices")) then - call get_node_value(doc, "write_matices", temp_str) + call get_node_value(doc, "write_matrices", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & cmfd_write_matrices = .true. From 58723fa0ecf1a90665032392c88b0ed556006e0f Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 10:11:36 -0400 Subject: [PATCH 16/22] write matrices when option is set --- src/cmfd_solver.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 49efe988dd..11989ff4f6 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -104,9 +104,7 @@ contains use error, only: fatal_error use global, only: cmfd, message, cmfd_shift, keff, & cmfd_ktol, cmfd_stol -#ifdef PETSC use global, only: cmfd_write_matrices -#endif logical :: adjoint @@ -150,8 +148,10 @@ contains ! Setup petsc for everything call loss % assemble() call prod % assemble() - call loss % write('loss.dat') - call prod % write('prod.dat') + if (cmfd_write_matrices) then + call loss % write('loss.dat') + call prod % write('prod.dat') + end if ! Set norms to 0 norm_n = ZERO From d1050fa194dcfeea2338837a5659c2ae565d1ba0 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 10 Sep 2014 09:30:46 -0400 Subject: [PATCH 17/22] indented 5 spaces for line continuation --- src/cmfd_input.F90 | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 61d6f6eaba..0a27b6a19c 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -153,7 +153,7 @@ contains call get_node_value(doc, "feedback", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_feedback = .true. + cmfd_feedback = .true. end if ! Set downscatter logical @@ -161,31 +161,31 @@ contains call get_node_value(doc, "downscatter", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_downscatter = .true. + cmfd_downscatter = .true. end if ! Set the solver type if (check_for_node(doc, "solver")) & - call get_node_value(doc, "solver", cmfd_solver_type) + call get_node_value(doc, "solver", cmfd_solver_type) ! Set monitoring if (check_for_node(doc, "snes_monitor")) then call get_node_value(doc, "snes_monitor", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_snes_monitor = .true. + cmfd_snes_monitor = .true. end if if (check_for_node(doc, "ksp_monitor")) then call get_node_value(doc, "ksp_monitor", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_ksp_monitor = .true. + cmfd_ksp_monitor = .true. end if if (check_for_node(doc, "power_monitor")) then call get_node_value(doc, "power_monitor", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_power_monitor = .true. + cmfd_power_monitor = .true. end if ! Output logicals @@ -193,7 +193,7 @@ contains call get_node_value(doc, "write_matrices", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_write_matrices = .true. + cmfd_write_matrices = .true. end if ! Run an adjoint calc @@ -201,36 +201,36 @@ contains call get_node_value(doc, "run_adjoint", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_run_adjoint = .true. + cmfd_run_adjoint = .true. end if ! Batch to begin cmfd if (check_for_node(doc, "begin")) & - call get_node_value(doc, "begin", cmfd_begin) + call get_node_value(doc, "begin", cmfd_begin) ! Tally during inactive batches if (check_for_node(doc, "inactive")) then call get_node_value(doc, "inactive", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & - cmfd_tally_on = .false. + cmfd_tally_on = .false. end if ! Inactive batch flush window if (check_for_node(doc, "inactive_flush")) & - call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1)) + call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1)) if (check_for_node(doc, "num_flushes")) & - call get_node_value(doc, "num_flushes", cmfd_inact_flush(2)) + call get_node_value(doc, "num_flushes", cmfd_inact_flush(2)) ! Last flush before active batches if (check_for_node(doc, "active_flush")) & - call get_node_value(doc, "active_flush", cmfd_act_flush) + call get_node_value(doc, "active_flush", cmfd_act_flush) ! Get display if (check_for_node(doc, "display")) & - call get_node_value(doc, "display", cmfd_display) + call get_node_value(doc, "display", cmfd_display) if (trim(cmfd_display) == 'dominance' .and. & - trim(cmfd_solver_type) /= 'power') then + trim(cmfd_solver_type) /= 'power') then message = 'Dominance Ratio only aviable with power iteration solver' call warning() cmfd_display = '' @@ -238,17 +238,17 @@ contains ! Read in spectral radius estimate and tolerances if (check_for_node(doc, "spectral")) & - call get_node_value(doc, "spectral", cmfd_spectral) + call get_node_value(doc, "spectral", cmfd_spectral) if (check_for_node(doc, "shift")) & - call get_node_value(doc, "shift", cmfd_shift) + call get_node_value(doc, "shift", cmfd_shift) if (check_for_node(doc, "ktol")) & - call get_node_value(doc, "ktol", cmfd_ktol) + call get_node_value(doc, "ktol", cmfd_ktol) if (check_for_node(doc, "stol")) & - call get_node_value(doc, "stol", cmfd_stol) + call get_node_value(doc, "stol", cmfd_stol) if (check_for_node(doc, "atoli")) & - call get_node_value(doc, "atoli", cmfd_atoli) + call get_node_value(doc, "atoli", cmfd_atoli) if (check_for_node(doc, "rtoli")) & - call get_node_value(doc, "rtoli", cmfd_rtoli) + call get_node_value(doc, "rtoli", cmfd_rtoli) ! Create tally objects call create_cmfd_tally(doc) @@ -421,7 +421,7 @@ contains call get_node_value(doc, "reset", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - t % reset = .true. + t % reset = .true. end if ! Set up mesh filter From 736f236a284d6d2cda9f6318b6a80dfeac85467d Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 10 Sep 2014 09:36:46 -0400 Subject: [PATCH 18/22] remove commented out code --- src/cmfd_solver.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 11989ff4f6..9a2ebceb1b 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -429,7 +429,6 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) its = igs -!print *, igs, err, w if (err < tol) exit ! Calculation new overrelaxation parameter @@ -569,8 +568,6 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) its = igs -!print *, igs, err, w -!read * if (err < tol) exit ! Calculation new overrelaxation parameter From 1c797abc38af1a774b182fd78e855ea03240f0a3 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 12 Sep 2014 09:30:05 -0400 Subject: [PATCH 19/22] added rest of intents and moved convergence logical --- src/cmfd_solver.F90 | 74 ++++++++++++++++++++++++------------------- src/matrix_header.F90 | 18 +++++------ src/vector_header.F90 | 4 +-- 3 files changed, 52 insertions(+), 44 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 9a2ebceb1b..69a2829039 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -11,7 +11,6 @@ module cmfd_solver private public :: cmfd_solver_execute - logical :: iconv ! did the problem converged real(8) :: k_n ! new k-eigenvalue real(8) :: k_o ! old k-eigenvalue real(8) :: k_s ! shift of eigenvalue @@ -38,11 +37,11 @@ module cmfd_solver subroutine linsolve(A, b, x, tol, i) import :: Matrix import :: Vector - type(Matrix) :: A - type(Vector) :: b - type(Vector) :: x - real(8) :: tol - integer :: i + type(Matrix), intent(inout) :: A + type(Vector), intent(inout) :: b + type(Vector), intent(inout) :: x + real(8), intent(in) :: tol + integer, intent(out) :: i end subroutine linsolve end interface @@ -56,7 +55,7 @@ contains use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve - logical, optional :: adjoint ! adjoint calc + logical, optional, intent(in) :: adjoint ! adjoint calc logical :: physical_adjoint = .false. @@ -106,7 +105,7 @@ contains cmfd_ktol, cmfd_stol use global, only: cmfd_write_matrices - logical :: adjoint + logical, intent(in) :: adjoint integer :: n ! problem size real(8) :: guess ! initial guess @@ -145,7 +144,7 @@ contains ! Fill in production matrix call build_prod_matrix(prod, adjoint=adjoint) - ! Setup petsc for everything + ! Finalize setup of CSR matrices call loss % assemble() call prod % assemble() if (cmfd_write_matrices) then @@ -209,11 +208,13 @@ contains subroutine execute_power_iter() use constants, only: ONE - use global, only: cmfd_atoli, cmfd_rtoli + use error, only: fatal_error + use global, only: cmfd_atoli, cmfd_rtoli, message integer :: i ! iteration counter integer :: innerits ! # of inner iterations integer :: totalits ! total number of inners + logical :: iconv ! did the problem converged real(8) :: atoli ! absolute minimum tolerance real(8) :: rtoli ! relative tolerance based on source conv real(8) :: toli ! the current tolerance of inners @@ -233,6 +234,12 @@ contains ! Begin power iteration do i = 1, 10000 + ! Check if reached iteration 10000 + if (i == 10000) then + message = 'Reached maximum iterations in CMFD power iteration solver.' + call fatal_error() + end if + ! Compute source vector call prod % vector_multiply(phi_o, s_o) @@ -255,7 +262,7 @@ contains s_o % val = s_o % val * k_lo ! Check convergence - call convergence(i, innerits) + call convergence(i, innerits, iconv) totalits = totalits + innerits ! Break loop if converged @@ -304,14 +311,15 @@ contains ! CONVERGENCE checks the convergence of the CMFD problem !=============================================================================== - subroutine convergence(iter, innerits) + subroutine convergence(iter, innerits, iconv) use constants, only: ONE, TINY_BIT use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV - integer :: iter ! outer iteration number - integer :: innerits ! inner iteration nubmer + integer, intent(in) :: iter ! outer iteration number + integer, intent(in) :: innerits ! inner iteration nubmer + logical, intent(out) :: iconv ! convergence logical ! Reset convergence flag iconv = .false. @@ -349,11 +357,11 @@ contains use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral - type(Matrix) :: A ! coefficient matrix - type(Vector) :: b ! right hand side vector - type(Vector) :: x ! unknown vector - real(8) :: tol ! tolerance on final error - integer :: its ! number of inner iterations + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations integer :: g ! group index integer :: i ! loop counter for x @@ -449,11 +457,11 @@ contains use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral - type(Matrix) :: A ! coefficient matrix - type(Vector) :: b ! right hand side vector - type(Vector) :: x ! unknown vector - real(8) :: tol ! tolerance on final error - integer :: its ! number of inner iterations + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations integer :: g ! group index integer :: i ! loop counter for x @@ -646,15 +654,15 @@ contains use global, only: cmfd, cmfd_coremap - integer :: i ! iteration counter for x - integer :: j ! iteration counter for y - integer :: k ! iteration counter for z - integer :: g ! iteration counter for groups - integer :: irow ! iteration counter over row (0 reference) - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups + integer, intent(out) :: i ! iteration counter for x + integer, intent(out) :: j ! iteration counter for y + integer, intent(out) :: k ! iteration counter for z + integer, intent(out) :: g ! iteration counter for groups + integer, intent(in) :: irow ! iteration counter over row (0 reference) + integer, intent(in) :: nx ! maximum number of x cells + integer, intent(in) :: ny ! maximum number of y cells + integer, intent(in) :: nz ! maximum number of z cells + integer, intent(in) :: ng ! maximum number of groups ! Check for core map if (cmfd_coremap) then diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 3f5aa204f7..3e3d4ad0db 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -368,11 +368,11 @@ contains subroutine matrix_search_indices(self, row, col, idx, found) - class(Matrix) :: self - integer :: row - integer :: col - integer :: idx - logical :: found + class(Matrix), intent(inout) :: self + integer, intent(in) :: row + integer, intent(in) :: col + integer, intent(out) :: idx + logical, intent(out) :: found integer :: j @@ -396,8 +396,8 @@ contains subroutine matrix_write(self, filename) - character(*) :: filename - class(Matrix) :: self + character(*), intent(in) :: filename + class(Matrix), intent(inout) :: self integer :: unit_ integer :: i @@ -421,8 +421,8 @@ contains subroutine matrix_copy(self, mattocopy) - class(Matrix) :: self - type(Matrix) :: mattocopy + class(Matrix), intent(inout) :: self + type(Matrix), intent(in) :: mattocopy ! Set n and nnz self % n_count = mattocopy % n_count diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 4d8cb2cf5c..755d374628 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -132,8 +132,8 @@ contains subroutine vector_copy(self, vectocopy) - class(Vector), target :: self - type(Vector) :: vectocopy + class(Vector), target, intent(inout) :: self + type(Vector), intent(in) :: vectocopy ! Preallocate vector if (.not.allocated(self % data)) allocate(self % data(vectocopy % n)) From 70785e62bddd2bee02744c085f7db1eb104d592f Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 12 Sep 2014 10:01:56 -0400 Subject: [PATCH 20/22] CMFD tests no longer removed if PETSc is not enabled --- src/CMakeLists.txt | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fc6c6630cc..acf7314bbc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -147,9 +147,7 @@ endif() # PETSc for CMFD functionality #=============================================================================== -set (PETSC_ENABLED FALSE) if(petsc) - set(PETSC_ENABLED TRUE) find_package(PETSc REQUIRED HINTS $ENV{PETSC_DIR}/conf) find_library(libpetsc petsc $ENV{PETSC_DIR}/lib) @@ -275,14 +273,6 @@ include(CTest) # Get a list of all the tests to run file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py) -# Check to see if PETSC is compiled for CMFD tests -if (NOT ${PETSC_ENABLED}) - file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd*.py) - foreach(cmfd_test in ${CMFD_TESTS}) - list(REMOVE_ITEM TESTS ${cmfd_test}) - endforeach(cmfd_test) -endif(NOT ${PETSC_ENABLED}) - # Check for MEM_CHECK and COVERAGE variables if (DEFINED ENV{MEM_CHECK}) set(MEM_CHECK $ENV{MEM_CHECK}) From fe3b178cd41e4120ec45031c6bd90be5abae2aa0 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 12 Sep 2014 10:35:11 -0400 Subject: [PATCH 21/22] allow power iteration CMFD to run without PETSc --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fc6c6630cc..8a9424e79d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -277,7 +277,7 @@ file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py) # Check to see if PETSC is compiled for CMFD tests if (NOT ${PETSC_ENABLED}) - file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd*.py) + file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd_jfnk.py) foreach(cmfd_test in ${CMFD_TESTS}) list(REMOVE_ITEM TESTS ${cmfd_test}) endforeach(cmfd_test) From c8e5f984ad29c1fddc67949101eddacd0d43b1bf Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 12 Sep 2014 13:28:26 -0400 Subject: [PATCH 22/22] re-added compiler directives around petsc error int --- src/vector_header.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 755d374628..05cc7e3cde 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -28,7 +28,9 @@ module vector_header # endif end type Vector +#ifdef PETSC integer :: petsc_err ! petsc error code +#endif contains