mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Merge branch 'develop' into case-insensitive
Conflicts: src/cmfd_input.F90
This commit is contained in:
Paul Romano
commit
192743fda4
10 changed files with 877 additions and 36 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
705
src/cmfd_solver.F90
Normal file
705
src/cmfd_solver.F90
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -12,5 +12,7 @@
|
|||
<display> dominance </display>
|
||||
<solver> power </solver>
|
||||
<feedback> true </feedback>
|
||||
<atoli> 1.e-15 </atoli>
|
||||
<rtoli> 1.e-20 </rtoli>
|
||||
|
||||
</cmfd>
|
||||
|
|
|
|||
|
|
@ -12,5 +12,7 @@
|
|||
<display> dominance </display>
|
||||
<solver> power </solver>
|
||||
<feedback> false </feedback>
|
||||
<atoli> 1.e-15 </atoli>
|
||||
<rtoli> 1.e-20 </rtoli>
|
||||
|
||||
</cmfd>
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue