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)
diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index 6234d9fcf9..35485a8945 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_input.F90 b/src/cmfd_input.F90
index 3be4eb66a6..5ab4caf480 100644
--- a/src/cmfd_input.F90
+++ b/src/cmfd_input.F90
@@ -153,7 +153,7 @@ contains
call get_node_value(doc, "feedback", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
- cmfd_feedback = .true.
+ cmfd_feedback = .true.
end if
! Set downscatter logical
@@ -161,39 +161,39 @@ contains
call get_node_value(doc, "downscatter", temp_str)
temp_str = to_lower(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)
temp_str = to_lower(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)
temp_str = to_lower(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)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
- cmfd_power_monitor = .true.
+ cmfd_power_monitor = .true.
end if
! 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)
temp_str = to_lower(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,41 +201,55 @@ contains
call get_node_value(doc, "run_adjoint", temp_str)
temp_str = to_lower(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)
temp_str = to_lower(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 = ''
end if
+ ! Read in spectral radius estimate and tolerances
+ if (check_for_node(doc, "spectral")) &
+ call get_node_value(doc, "spectral", cmfd_spectral)
+ if (check_for_node(doc, "shift")) &
+ call get_node_value(doc, "shift", cmfd_shift)
+ if (check_for_node(doc, "ktol")) &
+ call get_node_value(doc, "ktol", cmfd_ktol)
+ if (check_for_node(doc, "stol")) &
+ call get_node_value(doc, "stol", cmfd_stol)
+ if (check_for_node(doc, "atoli")) &
+ call get_node_value(doc, "atoli", cmfd_atoli)
+ if (check_for_node(doc, "rtoli")) &
+ call get_node_value(doc, "rtoli", cmfd_rtoli)
+
! Create tally objects
call create_cmfd_tally(doc)
@@ -407,7 +421,7 @@ contains
call get_node_value(doc, "reset", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
- t % reset = .true.
+ t % reset = .true.
end if
! Set up mesh filter
diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90
new file mode 100644
index 0000000000..69a2829039
--- /dev/null
+++ b/src/cmfd_solver.F90
@@ -0,0 +1,705 @@
+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
+
+ 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) :: 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
+ 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
+
+ ! CMFD linear solver interface
+ procedure(linsolve), pointer :: cmfd_linsolver => null()
+ abstract interface
+ subroutine linsolve(A, b, x, tol, i)
+ import :: Matrix
+ import :: Vector
+ 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
+
+contains
+
+!===============================================================================
+! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD
+!===============================================================================
+
+ subroutine cmfd_solver_execute(adjoint)
+
+ use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
+
+ logical, optional, intent(in) :: adjoint ! adjoint calc
+
+ logical :: physical_adjoint = .false.
+
+ ! 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 error, only: fatal_error
+ use global, only: cmfd, message, cmfd_shift, keff, &
+ cmfd_ktol, cmfd_stol
+ use global, only: cmfd_write_matrices
+
+ logical, intent(in) :: adjoint
+
+ integer :: n ! problem size
+ real(8) :: guess ! initial guess
+ real(8) :: dw ! eigenvalue shift
+
+ ! 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 = 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)
+
+ ! Fill in production matrix
+ call build_prod_matrix(prod, adjoint=adjoint)
+
+ ! Finalize setup of CSR matrices
+ call loss % assemble()
+ call prod % assemble()
+ if (cmfd_write_matrices) then
+ call loss % write('loss.dat')
+ call prod % write('prod.dat')
+ end if
+
+ ! Set norms to 0
+ 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
+
+ ! Set tolerances
+ ktol = cmfd_ktol
+ stol = cmfd_stol
+
+ end subroutine init_data
+
+!===============================================================================
+! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem
+!===============================================================================
+
+ subroutine compute_adjoint()
+
+ use error, only: fatal_error
+ use global, only: message
+#ifdef PETSC
+ 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
+#else
+ message = 'Adjoint calculations only allowed with PETSc'
+ call fatal_error()
+#endif
+
+ end subroutine compute_adjoint
+
+!===============================================================================
+! EXECUTE_POWER_ITER is the main power iteration routine
+! for the cmfd calculation
+!===============================================================================
+
+ subroutine execute_power_iter()
+
+ use constants, only: ONE
+ 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
+
+ ! Reset convergence flag
+ iconv = .false.
+
+ ! Set up tolerances
+ atoli = cmfd_atoli
+ rtoli = cmfd_rtoli
+ toli = rtoli*100._8
+
+ ! Perform shift
+ call wielandt_shift()
+ totalits = 0
+
+ ! 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)
+
+ ! Normalize source vector
+ s_o % val = s_o % val / k_lo
+
+ ! Compute new flux vector
+ call cmfd_linsolver(loss, s_o, phi_n, toli, innerits)
+
+ ! Compute new source vector
+ call prod % vector_multiply(phi_n, s_n)
+
+ ! 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_lo
+
+ ! Check convergence
+ call convergence(i, innerits, iconv)
+ totalits = totalits + innerits
+
+ ! Break loop if converged
+ if (iconv) exit
+
+ ! Record old values
+ phi_o % val = phi_n % val
+ k_o = k_n
+ k_lo = k_ln
+ norm_o = norm_n
+
+ ! Get new tolerance for inners
+ toli = max(atoli, rtoli*serr)
+
+ end do
+
+ end subroutine execute_power_iter
+
+!===============================================================================
+! WIELANDT SHIFT
+!===============================================================================
+
+ subroutine wielandt_shift()
+
+ use constants, only: ONE
+
+ integer :: irow ! row counter
+ integer :: icol ! col counter
+ integer :: jcol ! current col index in prod matrix
+
+ ! perform subtraction
+ 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) .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
+ end do COLS
+ end do ROWS
+
+ end subroutine wielandt_shift
+
+!===============================================================================
+! CONVERGENCE checks the convergence of the CMFD problem
+!===============================================================================
+
+ subroutine convergence(iter, innerits, iconv)
+
+ use constants, only: ONE, TINY_BIT
+ use global, only: cmfd_power_monitor, master
+ use, intrinsic :: ISO_FORTRAN_ENV
+
+ 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.
+
+ ! 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,T80,I0)') iter, k_n, kerr, &
+ serr, innerits
+ end if
+
+ end subroutine convergence
+
+!===============================================================================
+! CMFD_LINSOLVER_1g solves the CMFD linear system
+!===============================================================================
+
+ subroutine cmfd_linsolver_1g(A, b, x, tol, its)
+
+ use constants, only: ONE, ZERO
+ use global, only: cmfd, cmfd_spectral
+
+ 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
+ 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) :: w ! overrelaxation parameter
+ type(Vector) :: tmpx ! temporary solution vector
+
+ ! Set overrelaxation parameter
+ w = ONE
+
+ ! Dimensions
+ ng = 1
+ nx = cmfd % indices(1)
+ ny = cmfd % indices(2)
+ nz = cmfd % indices(3)
+ n = A % 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 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 = A % get_row(irow), didx - 1
+ tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
+ end do
+ 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)/A % 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)
+ its = igs
+ 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(A, b, x, tol, its)
+
+ use constants, only: ONE, ZERO
+ use global, only: cmfd, cmfd_spectral
+
+ 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
+ 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 :: 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
+ 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) :: x1 ! new g1 value of x
+ real(8) :: x2 ! new g2 value of x
+ real(8) :: err ! error in convergence of solution
+ real(8) :: w ! overrelaxation parameter
+ type(Vector) :: tmpx ! temporary solution vector
+
+ ! Set tolerance and overrelaxation parameter
+ w = ONE
+
+ ! Dimensions
+ ng = 2
+ nx = cmfd % indices(1)
+ ny = cmfd % indices(2)
+ nz = cmfd % indices(3)
+ n = A % 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, 2
+
+ ! 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 A % search_indices(irow, irow, d1idx, found)
+ call A % search_indices(irow + 1, irow + 1, d2idx, found)
+
+ ! Get block 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
+ 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 = A % get_row(irow), d1idx - 1
+ tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
+ end do
+ 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, A % get_row(irow + 1) - 1
+ tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
+ end do
+ 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
+ tmp1 = b % val(irow) - tmp1
+ tmp2 = b % val(irow + 1) - tmp2
+
+ ! Solve for new x
+ x1 = d11*tmp1 + d12*tmp2
+ x2 = d21*tmp1 + d22*tmp2
+
+ ! Perform overrelaxation
+ x % val(irow) = (ONE - w)*x % val(irow) + w*x1
+ x % val(irow + 1) = (ONE - w)*x % val(irow + 1) + w*x2
+
+ end do ROWS
+
+ end do REDBLACK
+
+ ! Check convergence
+ err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
+ its = igs
+ 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_2g
+
+!===============================================================================
+! 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
+#ifdef PETSC
+ call phi_n % write_petsc_binary(filename)
+#endif
+ end if
+
+ end subroutine extract_results
+
+!===============================================================================
+! 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, 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
+
+ ! 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
+!===============================================================================
+
+ subroutine finalize()
+
+ ! Destroy all objects
+ 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/global.F90 b/src/global.F90
index b136c2c166..0a780aeb62 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -371,6 +371,14 @@ module global
! CMFD display info
character(len=25) :: cmfd_display = 'balance'
+ ! Estimate of spectral radius of CMFD matrices and tolerances
+ real(8) :: cmfd_spectral = ZERO
+ 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
+ real(8) :: cmfd_rtoli = 1.e-5_8
+
! Information about state points to be written
integer :: n_state_points = 0
type(SetInt) :: statepoint_batch
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 1e988bfbae..d21c40edc1 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -780,12 +780,6 @@ contains
temp_str = to_lower(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
diff --git a/src/matrix_header.F90 b/src/matrix_header.F90
index 24307b10bb..3e3d4ad0db 100644
--- a/src/matrix_header.F90
+++ b/src/matrix_header.F90
@@ -20,19 +20,22 @@ 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 :: vector_multiply => matrix_vector_multiply
-#ifdef PETSC
- procedure :: transpose => matrix_transpose
+ 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
-#endif
+ procedure :: transpose => matrix_transpose
+# endif
end type matrix
#ifdef PETSC
@@ -359,4 +362,87 @@ 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), intent(inout) :: self
+ integer, intent(in) :: row
+ integer, intent(in) :: col
+ integer, intent(out) :: idx
+ logical, intent(out) :: 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(*), intent(in) :: filename
+ class(Matrix), intent(inout) :: 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
+
+!===============================================================================
+! MATRIX_COPY copies a matrix
+!===============================================================================
+
+ subroutine matrix_copy(self, mattocopy)
+
+ class(Matrix), intent(inout) :: self
+ type(Matrix), intent(in) :: 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/vector_header.F90 b/src/vector_header.F90
index e050981c76..05cc7e3cde 100644
--- a/src/vector_header.F90
+++ b/src/vector_header.F90
@@ -21,10 +21,11 @@ module vector_header
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
-#ifdef PETSC
+ procedure :: copy => vector_copy
+# ifdef PETSC
procedure :: setup_petsc => vector_setup_petsc
procedure :: write_petsc_binary => vector_write_petsc_binary
-#endif
+# endif
end type Vector
#ifdef PETSC
@@ -127,4 +128,28 @@ contains
end subroutine vector_write_petsc_binary
#endif
+!===============================================================================
+! VECTOR_COPY allocates a separate vector and copies
+!===============================================================================
+
+ subroutine vector_copy(self, vectocopy)
+
+ class(Vector), target, intent(inout) :: self
+ type(Vector), intent(in) :: 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
diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml
index 51def74088..390e558639 100644
--- a/tests/test_cmfd_feed/cmfd.xml
+++ b/tests/test_cmfd_feed/cmfd.xml
@@ -12,5 +12,7 @@
dominance
power
true
+ 1.e-15
+ 1.e-20
diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml
index c36a396043..4a4bf30108 100644
--- a/tests/test_cmfd_nofeed/cmfd.xml
+++ b/tests/test_cmfd_nofeed/cmfd.xml
@@ -12,5 +12,7 @@
dominance
power
false
+ 1.e-15
+ 1.e-20