crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

added intent attribute on all procedure args

Browse source
This commit is contained in:
Bryan Herman 2013-09-29 15:07:42 -04:00
parent 0fcde07bc5
commit 4907ed2df2
11 changed files with 221 additions and 220 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/cmfd_data.F90
View file

@ -767,15 +767,15 @@ contains
use global, only: cmfd, cmfd_hold_weights
real(8) :: get_reflector_albedo ! reflector albedo
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 :: l ! iteration counter for leakages
integer, intent(in) :: i ! iteration counter for x
integer, intent(in) :: j ! iteration counter for y
integer, intent(in) :: k ! iteration counter for z
integer, intent(in) :: g ! iteration counter for groups
integer, intent(in) :: l ! iteration counter for leakages
integer :: shift_idx ! parameter to shift index by +1 or -1
real(8) :: current(12) ! partial currents for all faces of mesh cell
real(8) :: albedo ! the albedo
integer :: shift_idx ! parameter to shift index by +1 or -1
real(8) :: current(12) ! partial currents for all faces of mesh cell
real(8) :: albedo ! the albedo
! Get partial currents from object
current = cmfd%current(:,g,i,j,k)
@ -798,7 +798,7 @@ contains
end function get_reflector_albedo
!===============================================================================
! FIX_NEUTRON_BALANCE
! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance
!===============================================================================
subroutine fix_neutron_balance()

41
src/cmfd_execute.F90
View file

@ -254,7 +254,7 @@ contains
end subroutine calc_fission_source
!===============================================================================
! CMFD_REWEIGHT
! CMFD_REWEIGHT performs weighting of particles in the source bank
!===============================================================================
subroutine cmfd_reweight(new_weights)
@ -271,17 +271,18 @@ contains
use mpi
#endif
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: i ! iteration counter
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: n_groups ! number of energy groups
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
logical :: new_weights ! calcualte new weights
logical, intent(in) :: new_weights ! calcualte new weights
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: i ! iteration counter
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: n_groups ! number of energy groups
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
type(StructuredMesh), pointer :: m ! point to mesh
real(8), allocatable :: egrid(:) ! energy grid
@ -389,13 +390,13 @@ contains
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 cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: ng ! maximum number of energy groups
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of cells in x direction
integer, intent(in) :: ny ! maximum number of cells in y direction
integer, intent(in) :: ng ! maximum number of energy groups
! Check if coremap is used
if (cmfd_coremap) then
@ -413,7 +414,7 @@ contains
end function get_matrix_idx
!===============================================================================
! CMFD_TALLY_RESET
! CMFD_TALLY_RESET resets all cmfd tallies
!===============================================================================
subroutine cmfd_tally_reset()

6
src/cmfd_header.F90
View file

@ -93,8 +93,8 @@ contains
subroutine allocate_cmfd(this, n_batches)
integer :: n_batches
type(cmfd_type) :: this
integer, intent(in) :: n_batches ! number of batches in calc
type(cmfd_type), intent(inout) :: this ! cmfd instance
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
@ -169,7 +169,7 @@ contains
subroutine deallocate_cmfd(this)
type(cmfd_type) :: this
type(cmfd_type), intent(inout) :: this ! cmfd instance
if (allocated(this % egrid)) deallocate(this % egrid)
if (allocated(this % totalxs)) deallocate(this % totalxs)

2
src/cmfd_input.F90
View file

@ -228,7 +228,7 @@ contains
integer :: ng ! number of energy groups (default 1)
integer :: n_filters ! number of filters
integer :: i_filter_mesh ! index for mesh filter
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: filename ! name of cmfd file
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(TallyFilter) :: filters(N_FILTER_TYPES) ! temporary filters

16
src/cmfd_jfnk_solver.F90
View file

@ -30,7 +30,7 @@ contains
use global, only: time_cmfdbuild, time_cmfdsolve
logical, optional :: adjoint ! adjoint calculation
logical, intent(in), optional :: adjoint ! adjoint calculation
! Check for adjoint
adjoint_calc = .false.
@ -84,8 +84,8 @@ contains
use constants, only: ZERO, ONE
use global, only: cmfd, cmfd_adjoint_type, current_batch
logical :: physical_adjoint
integer :: n
logical :: physical_adjoint ! physical adjoing calculation logical
integer :: n ! size of matrices
! Set up all matrices
call init_loss_matrix(loss)
@ -153,8 +153,8 @@ contains
subroutine init_jacobian_matrix()
integer :: nnz
integer :: n
integer :: nnz ! number of nonzeros in matrix
integer :: n ! dimension of matrix
! Get length of matrix and number of nonzeros total in loss matrix
nnz = loss % nnz
@ -179,7 +179,7 @@ contains
use constants, only: ONE
type(Vector) :: x
type(Vector), intent(in) :: x ! solution vector
integer :: i ! loop counter for jacobian rows
integer :: jjac ! loop counter for jacobian cols
@ -267,8 +267,8 @@ contains
use global, only: cmfd_write_matrices
type(Vector) :: x
type(Vector) :: res
type(Vector), intent(in) :: x ! solution vector
type(Vector), intent(out) :: res ! residual vector
character(len=25) :: filename
integer :: n

112
src/cmfd_loss_operator.F90
View file

@ -16,7 +16,7 @@ contains
subroutine init_loss_matrix(loss_matrix)
type(Matrix) :: loss_matrix ! cmfd loss matrix
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
@ -68,12 +68,12 @@ contains
function preallocate_loss_matrix(nx, ny, nz, ng, n) result(nnz)
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
integer, intent(in) :: n ! total length of matrix
integer :: nnz ! number of nonzeros
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
integer, intent(in) :: n ! total length of matrix
integer :: nnz ! number of nonzeros
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
@ -181,39 +181,39 @@ contains
subroutine build_loss_matrix(loss_matrix, adjoint)
type(Matrix) :: loss_matrix ! cmfd loss matrix
logical, optional :: adjoint ! set up the adjoint
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
logical, intent(in), optional :: adjoint ! set up the adjoint
integer :: nxyz(3,2) ! single vector containing bound. locations
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 :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
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 :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint calculation?
real(8) :: totxs ! total macro cross section
real(8) :: scattxsgg ! scattering macro cross section g-->g
real(8) :: scattxshg ! scattering macro cross section h-->g
real(8) :: dtilde(6) ! finite difference coupling parameter
real(8) :: dhat(6) ! nonlinear coupling parameter
real(8) :: hxyz(3) ! cell lengths in each direction
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
integer :: nxyz(3,2) ! single vector containing bound. locations
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 :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
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 :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint calculation?
real(8) :: totxs ! total macro cross section
real(8) :: scattxsgg ! scattering macro cross section g-->g
real(8) :: scattxshg ! scattering macro cross section h-->g
real(8) :: dtilde(6) ! finite difference coupling parameter
real(8) :: dhat(6) ! nonlinear coupling parameter
real(8) :: hxyz(3) ! cell lengths in each direction
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
! Check for adjoint
adjoint_calc = .false.
@ -366,14 +366,14 @@ contains
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
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
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
! Check if coremap is used
if (cmfd_coremap) then
@ -396,15 +396,15 @@ contains
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
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

12
src/cmfd_power_solver.F90
View file

@ -41,11 +41,11 @@ contains
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
real(8), intent(in), optional :: k_tol ! tolerance on keff
real(8), intent(in), optional :: s_tol ! tolerance on source
logical, intent(in), optional :: adjoint ! adjoint calc
logical :: physical_adjoint = .false.
logical :: physical_adjoint = .false. ! physical adjoint default false
! Set tolerances if present
if (present(k_tol)) ktol = k_tol
@ -104,7 +104,7 @@ contains
use constants, only: ONE, ZERO
use global, only: cmfd_write_matrices
logical :: adjoint
logical, intent(in) :: adjoint ! adjoint calcualtion
integer :: n ! problem size
real(8) :: guess ! initial guess
@ -238,7 +238,7 @@ contains
use global, only: cmfd_power_monitor, master
use, intrinsic :: ISO_FORTRAN_ENV
integer :: iter ! iteration number
integer, intent(in) :: iter ! iteration number
! Reset convergence flag
iconv = .false.

40
src/cmfd_prod_operator.F90
View file

@ -16,7 +16,7 @@ contains
subroutine init_prod_matrix(prod_matrix)
type(Matrix) :: prod_matrix
type(Matrix), intent(inout) :: prod_matrix ! production matrix
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
@ -50,8 +50,8 @@ contains
subroutine build_prod_matrix(prod_matrix, adjoint)
type(Matrix) :: prod_matrix
logical, optional :: adjoint
type(Matrix), intent(inout) :: prod_matrix ! production matrix
logical, intent(in), optional :: adjoint ! adjoint calculation logical
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
@ -135,14 +135,14 @@ contains
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
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
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
! check if coremap is used
if (cmfd_coremap) then
@ -165,15 +165,15 @@ contains
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
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

86
src/matrix_header.F90
View file

@ -43,9 +43,9 @@ contains
subroutine matrix_create(self, n, nnz)
integer :: n
integer :: nnz
class(Matrix) :: self
integer, intent(in) :: n ! dimension of matrix
integer, intent(in) :: nnz ! number of nonzeros
class(Matrix), intent(inout) :: self ! matrix instance
! Preallocate vectors
if (.not.allocated(self % row)) allocate(self % row(n+1))
@ -71,7 +71,7 @@ contains
subroutine matrix_destroy(self)
class(Matrix) :: self
class(Matrix), intent(inout) :: self ! matrix instance
#ifdef PETSC
if (self % petsc_active) call MatDestroy(self % petsc_mat, petsc_err)
@ -89,9 +89,9 @@ contains
subroutine matrix_add_value(self, col, val)
integer :: col
real(8) :: val
class(Matrix) :: self
integer, intent(in) :: col ! col location in matrix
real(8), intent(in) :: val ! value to store in matrix
class(Matrix), intent(inout) :: self ! matrix instance
! Record the data
self % col(self % nz_count) = col
@ -112,7 +112,7 @@ contains
subroutine matrix_new_row(self)
class(Matrix) :: self
class(Matrix), intent(inout) :: self ! matrix instance
! Record the current number of nonzeros
self % row(self % n_count) = self % nz_count
@ -132,7 +132,7 @@ contains
subroutine matrix_assemble(self)
class(Matrix) :: self
class(Matrix), intent(inout) :: self ! matrix instance
integer :: i
integer :: first
@ -158,12 +158,12 @@ contains
recursive subroutine sort_csr(col, val, first, last)
integer :: col(:)
integer :: first
integer :: last
real(8) :: val(:)
integer :: col(:) ! column vector to sort
integer :: first ! first value in sort
integer :: last ! last value in sort
real(8) :: val(:) ! value vector to be sorted like columns
integer :: mid
integer :: mid ! midpoint value
if (first < last) then
call split(col, val, first, last, mid)
@ -179,18 +179,18 @@ contains
subroutine split(col, val, low, high, mid)
integer :: col(:)
integer :: low
integer :: high
integer :: mid
real(8) :: val(:)
integer :: col(:) ! column vector to sort
integer :: low ! low index to sort
integer :: high ! high index to sort
integer :: mid ! middle of sort
real(8) :: val(:) ! value vector to be sorted like columns
integer :: left
integer :: right
integer :: iswap
integer :: pivot
real(8) :: rswap
real(8) :: val0
integer :: left ! contains left value in sort
integer :: right ! contains right value in sort
integer :: iswap ! temporary interger swap
integer :: pivot ! pivotting variable for columns
real(8) :: rswap ! temporary real swap
real(8) :: val0 ! pivot for value vector
left = low
right = high
@ -232,14 +232,14 @@ contains
end subroutine split
!===============================================================================
! MATRIX_GET_ROW
! MATRIX_GET_ROW is a method to get row and checks for PETSc C indexing use
!===============================================================================
function matrix_get_row(self, i) result(row)
class(Matrix) :: self
integer :: i
integer :: row
class(Matrix), intent(in) :: self ! matrix instance
integer, intent(in) :: i ! row to get
integer :: row ! index in col where row begins
row = self % row(i)
@ -248,14 +248,14 @@ contains
end function matrix_get_row
!===============================================================================
! MATRIX_GET_COL
! MATRIX_GET_COL is a method to get column and checks for PETSc C index use
!===============================================================================
function matrix_get_col(self, i) result(col)
class(Matrix) :: self
integer :: i
integer :: col
class(Matrix), intent(in) :: self ! matrix instance
integer, intent(in) :: i ! index from row vector
integer :: col ! the actual column
col = self % col(i)
@ -269,7 +269,7 @@ contains
subroutine matrix_setup_petsc(self)
class(Matrix) :: self
class(Matrix), intent(inout) :: self ! matrix instance
! change indices to c notation
self % row = self % row - 1
@ -292,11 +292,11 @@ contains
subroutine matrix_write_petsc_binary(self, filename)
character(*) :: filename
class(Matrix) :: self
character(*), intent(in) :: filename ! file name to write to
class(Matrix), intent(in) :: self ! matrix instance
#ifdef PETSC
type(PetscViewer) :: viewer
type(PetscViewer) :: viewer ! a petsc viewer instance
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
FILE_MODE_WRITE, viewer, petsc_err)
@ -312,7 +312,7 @@ contains
subroutine matrix_transpose(self)
class(Matrix) :: self
class(Matrix), intent(inout) :: self ! matrix instance
#ifdef PETSC
call MatTranspose(self % petsc_mat, MAT_REUSE_MATRIX, self % petsc_mat, &
@ -330,12 +330,12 @@ contains
use constants, only: ZERO
use vector_header, only: Vector
class(Matrix) :: self
type(Vector) :: vec_in
type(Vector) :: vec_out
class(Matrix), intent(in) :: self ! matrix instance
type(Vector), intent(in) :: vec_in ! vector to multiply matrix against
type(Vector), intent(inout) :: vec_out ! resulting vector
integer :: i
integer :: j
integer :: i ! row iteration counter
integer :: j ! column iteration counter
! Begin loop around rows
ROWS: do i = 1, self % n

82
src/solver_interface.F90
View file

@ -16,8 +16,8 @@ module solver_interface
! GMRES solver type
type, public :: GMRESSolver
#ifdef PETSC
type(ksp) :: ksp_
type(pc) :: pc_
type(ksp) :: ksp_ ! Krylov linear solver instance
type(pc) :: pc_ ! Preconditioner instance
#endif
contains
#ifdef PETSC
@ -37,11 +37,11 @@ module solver_interface
! JFNK solver type
type, public :: JFNKSolver
#ifdef PETSC
type(ksp) :: ksp_
type(pc) :: pc_
type(snes) :: snes_
type(mat) :: jac_mf
integer :: ls
type(ksp) :: ksp_ ! Krylov linear solver instance
type(pc) :: pc_ ! Preconditioner instance
type(snes) :: snes_ ! Nonlinear solver instance
type(mat) :: jac_mf ! Matrix free jacobian instance
integer :: ls ! Line search instance
#endif
contains
#ifdef PETSC
@ -56,13 +56,13 @@ module solver_interface
abstract interface
subroutine res_interface(x, res)
import :: Vector
type(Vector) :: x
type(Vector) :: res
type(Vector), intent(in) :: x ! solution vector
type(Vector), intent(out) :: res ! residual vector
end subroutine res_interface
subroutine jac_interface(x)
import :: Vector
type(Vector) :: x
type(Vector), intent(in) :: x ! solution vector
end subroutine jac_interface
end interface
@ -79,7 +79,7 @@ contains
subroutine petsc_gmres_create(self)
class(GMRESSolver) :: self
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
integer :: ilu_levels = 5
real(8) :: rtol = 1.0e-10_8
@ -102,9 +102,9 @@ contains
subroutine petsc_gmres_set_oper(self, prec_mat, mat_in)
class(GMRESSolver) :: self
type(Matrix) :: prec_mat
type(Matrix) :: mat_in
class(GMRESSolver), intent(inout) :: self ! GMRES solver instanace
type(Matrix), intent(inout) :: prec_mat ! preconditioner matrix
type(Matrix), intent(inout) :: mat_in ! coefficient matrix
call KSPSetOperators(self % ksp_, mat_in % petsc_mat, prec_mat % petsc_mat, &
SAME_NONZERO_PATTERN, petsc_err)
@ -118,7 +118,7 @@ contains
subroutine petsc_gmres_destroy(self)
class(GMRESSolver) :: self
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
call KSPDestroy(self % ksp_, petsc_err)
@ -130,9 +130,9 @@ contains
subroutine petsc_gmres_solve(self, b, x)
class(GMRESSolver) :: self
type(Vector) :: b
type(Vector) :: x
class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! solution vector
call KSPSolve(self % ksp_, b % petsc_vec, x % petsc_vec, petsc_err)
@ -144,7 +144,7 @@ contains
subroutine petsc_jfnk_create(self)
class(JFNKSolver) :: self
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
! Turn on mf_operator option for matrix free jacobian
call PetscOptionsSetValue("-snes_mf_operator", "TRUE", petsc_err)
@ -176,7 +176,7 @@ contains
subroutine petsc_jfnk_destroy(self)
class(JFNKSolver) :: self
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
call MatDestroy(self % jac_mf, petsc_err)
call SNESDestroy(self % snes_, petsc_err)
@ -189,10 +189,10 @@ contains
subroutine petsc_jfnk_set_functions(self, ctx, res, jac_prec)
class(JFNKSolver) :: self
type(Jfnk_ctx) :: ctx
type(Vector) :: res
type(Matrix) :: jac_prec
class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
type(Vector), intent(inout) :: res ! residual vector
type(Matrix), intent(inout) :: jac_prec ! preconditioner matrix
! Set residual procedure
call SNESSetFunction(self % snes_, res % petsc_vec, &
@ -217,8 +217,8 @@ contains
subroutine petsc_jfnk_solve(self, xvec)
class(JFNKSolver) :: self
type(Vector) :: xvec
class(JFNKSolver), intent(inout) :: self ! JFNK instance
type(Vector), intent(inout) :: xvec ! solution vector
call SNESSolve(self % snes_, PETSC_NULL_DOUBLE, xvec % petsc_vec, petsc_err)
@ -230,14 +230,14 @@ contains
subroutine petsc_jfnk_compute_residual(snes_, x, res, ctx, ierr)
type(snes) :: snes_
type(vec) :: x
type(vec) :: res
integer :: ierr
type(Jfnk_ctx) :: ctx
type(snes), intent(inout) :: snes_ ! PETSc SNES object
type(vec), intent(inout) :: x ! PETSc solution vector
type(vec), intent(inout) :: res ! PETSc residual vector
integer, intent(inout) :: ierr ! error code
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
type(Vector) :: xvec
type(Vector) :: resvec
type(Vector) :: xvec ! solution vector
type(Vector) :: resvec ! residual vector
! We need to use the x and res that come from PETSc because the pointer
! location changes and therefore we can not use module variables
@ -263,15 +263,15 @@ contains
subroutine petsc_jfnk_compute_jacobian(snes_, x, jac_mf, jac_prec, flag, &
ctx, ierr)
type(snes) :: snes_
type(vec) :: x
type(mat) :: jac_mf
type(mat) :: jac_prec
integer :: flag
type(Jfnk_ctx) :: ctx
integer :: ierr
type(snes), intent(inout) :: snes_ ! PETSc snes instance
type(vec), intent(inout) :: x ! PETSc solution vector
type(mat), intent(inout) :: jac_mf ! PETSc matrix free jacobian
type(mat), intent(inout) :: jac_prec ! PETSc matrix jacobian precond.
integer, intent(inout) :: flag ! unused madatory flag
type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
integer, intent(inout) :: ierr ! error code
type(Vector) :: xvec
type(Vector) :: xvec ! solution vector
! Again, we use the vector that comes from Petsc to build the Jacobian
! matrix

26
src/vector_header.F90
View file

@ -14,9 +14,9 @@ module vector_header
real(8), allocatable :: data(:) ! where vector data is stored
real(8), pointer :: val(:) ! pointer to vector data
# ifdef PETSC
type(vec) :: petsc_vec
type(vec) :: petsc_vec ! PETSc vector
# endif
logical :: petsc_active
logical :: petsc_active ! Logical if PETSc is being used
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
@ -25,7 +25,7 @@ module vector_header
procedure :: write_petsc_binary => vector_write_petsc_binary
end type Vector
integer :: petsc_err
integer :: petsc_err ! petsc error code
contains
@ -35,8 +35,8 @@ contains
subroutine vector_create(self, n)
integer :: n
class(Vector), target :: self
integer, intent(in) :: n ! size of vector
class(Vector), intent(inout), target :: self ! vector instance
! Preallocate vector
if (.not.allocated(self % data)) allocate(self % data(n))
@ -59,7 +59,7 @@ contains
subroutine vector_destroy(self)
class(Vector) :: self
class(Vector), intent(inout) :: self ! vector instance
#ifdef PETSC
if (self % petsc_active) call VecDestroy(self % petsc_vec, petsc_err)
@ -76,9 +76,9 @@ contains
subroutine vector_add_value(self, idx, val)
integer :: idx
real(8) :: val
class(Vector) :: self
integer, intent(in) :: idx ! index location in vector
real(8), intent(in) :: val ! value to add
class(Vector), intent(inout) :: self ! vector instance
self % val(idx) = val
@ -90,7 +90,7 @@ contains
subroutine vector_setup_petsc(self)
class(Vector) :: self
class(Vector), intent(inout) :: self ! vector instance
! Link to PETSc
#ifdef PETSC
@ -109,11 +109,11 @@ contains
subroutine vector_write_petsc_binary(self, filename)
character(*) :: filename
class(Vector) :: self
character(*), intent(in) :: filename ! name of file to write to
class(Vector), intent(in) :: self ! vector instance
#ifdef PETSC
type(PetscViewer) :: viewer
type(PetscViewer) :: viewer ! PETSc viewer instance
call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
FILE_MODE_WRITE, viewer, petsc_err)