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

created jacobian operator routine, all is parallelized except for building jacobian routine, iterfaces are updated between subroutines

Browse source
This commit is contained in:
Bryan Herman 2012-02-12 12:36:18 -08:00
parent 6ac8aec8c0
commit d62bbbf9a2
7 changed files with 402 additions and 178 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/DEPENDENCIES
View file

@ -34,6 +34,8 @@ cmfd_input.o: string.o
cmfd_input.o: tally_header.o
cmfd_input.o: xml-fortran/templates/cmfd_t.o
cmfd_jacobian_operator.o: global.o
cmfd_loss_operator.o: global.o
cmfd_output.o: global.o
@ -49,6 +51,7 @@ cmfd_prod_operator.o: global.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: cmfd_jacobian_operator.o
cmfd_snes_solver.o: cmfd_loss_operator.o
cmfd_snes_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: cmfd_slepc_solver.o

2
src/Makefile
View file

@ -18,7 +18,7 @@ include OBJECTS
F90X = /home/bherman/petsc/release/bin/mpif90
F90XFLAGS = -Wall -Wno-unused-variable
COMPILER = petsc
DEBUG = yes
DEBUG = no
PROFILE = yes
OPTIMIZE = no
USE_MPI = no

1
src/OBJECTS
View file

@ -6,6 +6,7 @@ cmfd_data.o \
cmfd_execute.o \
cmfd_header.o \
cmfd_input.o \
cmfd_jacobian_operator.o \
cmfd_loss_operator.o \
cmfd_output.o \
cmfd_power_solver.o \

4
src/cmfd_execute.F90
View file

@ -41,7 +41,7 @@ contains
if (master) then
! begin timer
call timer_start(time_cmfd)
! call timer_start(time_cmfd)
! set up cmfd
if(.not. cmfd_only) call set_up_cmfd()
@ -57,7 +57,7 @@ contains
#endif
! stop timer
call timer_stop(time_cmfd)
! call timer_stop(time_cmfd)
! write vtk file
!if(.not. cmfd_only) call write_cmfd_vtk()

338
src/cmfd_jacobian_operator.F90 Normal file
View file

@ -0,0 +1,338 @@
module cmfd_jacobian_operator
use cmfd_loss_operator, only: loss_operator,init_M_operator, &
& build_loss_matrix,destroy_M_operator
use cmfd_prod_operator, only: prod_operator,init_F_operator, &
& build_prod_matrix,destroy_F_operator
#ifdef PETSC
implicit none
private
public :: init_J_operator,build_jacobian_matrix,destroy_J_operator
#include <finclude/petsc.h90>
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 :: ierr ! petsc error code
type, public :: jacobian_operator
Mat :: J ! petsc matrix for neutronic prod operator
integer :: n ! dimensions of matrix
integer :: nnz ! max number of nonzeros
integer :: localn ! local size on proc
integer, allocatable :: d_nnz(:) ! vector of diagonal preallocation
integer, allocatable :: o_nnz(:) ! vector of off-diagonal preallocation
end type jacobian_operator
type, public :: operators
type(loss_operator) :: loss
type(prod_operator) :: prod
end type operators
contains
!==============================================================================
! INIT_J_OPERATOR
!==============================================================================
subroutine init_J_operator(this,ctx)
type(jacobian_operator) :: this
type(operators) :: ctx
! get indices
call get_J_indices(this)
! get preallocation
call preallocate_jacobian_matrix(this,ctx)
! set up M operator
call MatCreateMPIAIJ(PETSC_COMM_WORLD,this%localn,this%localn,PETSC_DECIDE,&
& PETSC_DECIDE,PETSC_NULL_INTEGER,this%d_nnz,PETSC_NULL_INTEGER,this%o_nnz, &
& this%J,ierr)
call MatSetOption(this%J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr)
call MatSetOption(this%J,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr)
end subroutine init_J_operator
!==============================================================================
! GET_J_INDICES
!==============================================================================
subroutine get_J_indices(this)
use global, only: cmfd,cmfd_coremap
type(jacobian_operator) :: this
! get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! get number of nonzeros
this%nnz = 7 + ng - 1
! calculate dimensions of matrix
if (cmfd_coremap) then
this%n = cmfd % mat_dim * ng
else
this%n = nx*ny*nz*ng
end if
! add 1 for eigenvalue row
this%n = this%n + 1
end subroutine get_J_indices
!===============================================================================
! PREALLOCATE_JACOBIAN_MATRIX
!===============================================================================
subroutine preallocate_jacobian_matrix(this,ctx)
use global, only: cmfd,cmfd_coremap
type(jacobian_operator) :: this
type(operators) :: ctx
integer :: rank ! rank of processor
integer :: sizen ! number of procs
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 :: n ! the extent of the matrix
integer :: irow ! row counter
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 :: nxyz(3,2) ! single vector containing bound. locations
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: row_start ! index of local starting row
integer :: row_end ! index of local final row
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
! get rank and max rank of procs
call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizen,ierr)
! get local problem size
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(sizen)
! add 1 more if less proc id is less than mod
if (rank < mod(n,sizen)) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n,sizen)) then
row_start = rank*(n/sizen+1)
else
row_start = min(mod(n,sizen)*(n/sizen+1)+(rank - mod(n,sizen))*(n/sizen),n)
end if
! determine local final row
row_end = row_start + this%localn - 1
! allocate counters
if (.not. allocated(this%d_nnz)) allocate(this%d_nnz(row_start:row_end))
if (.not. allocated(this%o_nnz)) allocate(this%o_nnz(row_start:row_end))
this % d_nnz = 0
this % o_nnz = 0
! start with pattern from loss matrix
this % d_nnz = ctx%loss%d_nnz
this % o_nnz = ctx%loss%o_nnz
! append -F*phi term for last processor will take care of 1 for lambda
if (row_end == n - 1) then
this%d_nnz = this%d_nnz + 1
else
this%o_nnz = this%o_nnz + 1
end if
! do last row which has all filled (already did lower left corner above)
if (row_end == n - 1) then
this % d_nnz(row_end) = this % d_nnz(row_end) + (row_end - row_start)
this % o_nnz(row_end) = this % o_nnz(row_end) + ((this%n-1) - (row_end - &
& row_start))
end if
end subroutine preallocate_jacobian_matrix
!===============================================================================
! BUILD_JACOBIAN_MATRIX creates the matrix representing loss of neutrons
!===============================================================================
subroutine build_jacobian_matrix(snes,x,jac,jac_prec,flag,ctx,ierr)
! formal variables
SNES :: snes ! the snes context
Vec :: x ! the solution vector
Mat :: jac ! the jacobian matrix
Mat :: jac_prec ! the jacobian preconditioner
MatStructure :: flag ! not used
type(operators) :: ctx ! not used
integer :: ierr ! petsc error flag
! local varibles
Vec :: phi ! flux vector
Vec :: source ! source vector
integer :: n ! problem size
integer :: k ! implied do loop counter
integer :: ncols ! number of nonzeros in cols
integer :: irow ! row counter
integer :: row_start! starting local row on process
integer :: row_end ! ending local row on process
integer, allocatable :: nnzv(:) ! vector of number of nonzeros for jac
integer, allocatable :: cols(:) ! vector of column numbers
real(8) :: lambda ! eigenvalue
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: sptr(:) ! pointer to source vector
real(8), allocatable :: vals(:) ! vector of row values
! create operators
call build_loss_matrix(ctx%loss)
call build_prod_matrix(ctx%prod)
! get problem size
n = ctx%loss%n
! get local size on each processor
call MatGetOwnershipRange(jac_prec,row_start,row_end,ierr)
! allocate cols and initialize to zero
if (.not. allocated(cols)) allocate(cols(maxval(ctx%loss%d_nnz+ctx%loss%o_nnz)))
if (.not. allocated(vals)) allocate(vals(maxval(ctx%loss%d_nnz+ctx%loss%o_nnz)))
cols = 0
vals = 0.0_8
! get pointers to residual vector
call VecGetArrayF90(x,xptr,ierr)
! create petsc vector for flux
call VecCreate(PETSC_COMM_SELF,phi,ierr)
call VecSetSizes(phi,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(phi,ierr)
! extract flux and eigenvalue
call VecPlaceArray(phi,xptr,ierr)
lambda = xptr(n+1)
! compute math (M-lambda*F) M is overwritten here
call MatAXPY(ctx%loss%M,-1.0_8*lambda,ctx%prod%F,SUBSET_NONZERO_PATTERN,ierr)
! create tmp petsc vector for source
call VecCreate(PETSC_COMM_SELF,source,ierr)
call VecSetSizes(source,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(source,ierr)
! perform math (-F*phi --> source)
call MatMult(ctx%prod%F,phi,source,ierr)
call VecScale(source,-1.0_8,ierr)
! get pointer to source
call VecGetArrayF90(source,sptr,ierr)
! begin loop to insert things into matrix
do irow = row_start,row_end - 1
! don't do last row
if (irow == n) cycle
! get row of matrix
call MatGetRow(ctx%loss%M,irow,ncols,cols,vals,ierr)
! set that row to Jacobian matrix
call MatSetValues(jac_prec,1,irow,ncols,cols(1:ncols),vals,INSERT_VALUES,ierr)
! restore the row
call MatRestoreRow(ctx%loss%M,irow,ncols,cols,vals,ierr)
! insert source value
call MatSetValue(jac_prec,irow,n,sptr(irow+1),INSERT_VALUES,ierr)
end do
! set values in last row of matrix
if (row_end - 1 == n) then
call MatSetValues(jac_prec,1,n,n,(/(k,k=0,n-1)/),xptr(1:n),INSERT_VALUES,ierr)
call MatSetValue(jac_prec,n,n,1.0_8,INSERT_VALUES,ierr)
end if
! assemble matrix
call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
! reset all vectors
call VecResetArray(phi,ierr)
! restore all vectors
call VecRestoreArrayF90(x,xptr,ierr)
call VecRestoreArrayF90(source,sptr,ierr)
! destroy all temporary objects
call VecDestroy(phi,ierr)
call VecDestroy(source,ierr)
! deallocate all temporary space
if (allocated(cols)) deallocate(cols)
if (allocated(vals)) deallocate(vals)
end subroutine build_jacobian_matrix
!===============================================================================
! PRINT_J_OPERATOR
!===============================================================================
subroutine print_J_operator(this)
use global, only: path_input
PetscViewer :: viewer
type(jacobian_operator) :: this
! write out matrix in binary file (debugging)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD,trim(path_input)//'jacobian.bin'&
& ,FILE_MODE_WRITE,viewer,ierr)
call MatView(this%J,viewer,ierr)
call PetscViewerDestroy(viewer,ierr)
end subroutine print_J_operator
!==============================================================================
! DESTROY_J_OPERATOR
!==============================================================================
subroutine destroy_J_operator(this)
type(jacobian_operator) :: this
! deallocate matrix
call MatDestroy(this%J,ierr)
! deallocate other parameters
if (allocated(this%d_nnz)) deallocate(this%d_nnz)
if (allocated(this%o_nnz)) deallocate(this%o_nnz)
end subroutine destroy_J_operator
#endif
end module cmfd_jacobian_operator

23
src/cmfd_slepc_solver.F90
View file

@ -31,16 +31,26 @@ contains
subroutine cmfd_slepc_execute()
use timing
use global, only:time_cmfd,master
call timer_start(time_cmfd)
! initialize data
call init_data()
! initialize solver
call init_solver()
call timer_stop(time_cmfd)
if(master) print *,'Init Time:',time_cmfd%elapsed
call timer_reset(time_cmfd)
call timer_start(time_cmfd)
! build operators
call build_loss_matrix(loss)
call build_prod_matrix(prod)
call timer_stop(time_cmfd)
if(master) print *,'Build Time',time_cmfd%elapsed
call timer_reset(time_cmfd)
call timer_start(time_cmfd)
! set operators to EPS object
call EPSSetOperators(eps,prod%F,loss%M,ierr)
@ -49,10 +59,15 @@ contains
! solve the system
call EPSSolve(eps,ierr)
call timer_stop(time_cmfd)
if(master) print *,'Solve Time:',time_cmfd%elapsed
call timer_reset(time_cmfd)
call timer_start(time_cmfd)
! extracts results to cmfd object
call extract_results()
call timer_stop(time_cmfd)
if(master) print *,'Extraction Time:',time_cmfd%elapsed
call timer_reset(time_cmfd)
! deallocate all slepc data
call finalize()

209
src/cmfd_snes_solver.F90
View file

@ -1,21 +1,24 @@
module cmfd_snes_solver
#ifdef PETSC
use cmfd_loss_operator, only: loss_operator,init_M_operator, &
& build_loss_matrix,destroy_M_operator
use cmfd_prod_operator, only: prod_operator,init_F_operator, &
& build_prod_matrix,destroy_F_operator
use cmfd_slepc_solver, only: cmfd_slepc_execute
use cmfd_loss_operator, only: loss_operator,init_M_operator, &
& build_loss_matrix,destroy_M_operator
use cmfd_prod_operator, only: prod_operator,init_F_operator, &
& build_prod_matrix,destroy_F_operator
use cmfd_jacobian_operator, only: jacobian_operator,init_J_operator, &
& build_jacobian_matrix,destroy_J_operator, &
& operators
use cmfd_slepc_solver, only: cmfd_slepc_execute
implicit none
#include <finclude/petsc.h90>
type(loss_operator) :: loss
type(prod_operator) :: prod
type(jacobian_operator) :: jac_prec
type(operators) :: ctx
Mat :: jac ! jacobian matrix
Mat :: jac_prec ! preconditioned jacobian matrix
! Mat :: jac_prec ! preconditioned jacobian matrix
Vec :: resvec ! residual vector
Vec :: xvec ! results
KSP :: ksp ! linear solver context
@ -70,12 +73,13 @@ print *,'solving system'
integer, allocatable :: nnzv(:) ! vector of number of nonzeros in jacobian
real(8), pointer :: xptr(:) ! solution pointer
! set up matrices
call init_M_operator(loss)
call init_F_operator(prod)
! set up operator matrices
call init_M_operator(ctx%loss)
call init_F_operator(ctx%prod)
call init_J_operator(jac_prec,ctx)
! get problem size
n = loss%n
n = jac_prec%n - 1
! create PETSc vectors
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n+1,resvec,ierr)
@ -87,7 +91,7 @@ print *,'solving system'
if (row_end == n + 1) row_end = n
! set flux in guess
call VecSetValues(xvec,row_end-row_start-1,(/(k,k=row_start,row_end-1)/), &
call VecSetValues(xvec,row_end-row_start-1,(/(k,k=row_start,row_end-1)/), &
& cmfd%phi(row_start+1:row_end),INSERT_VALUES,ierr)
call VecAssemblyBegin(xvec,ierr)
call VecAssemblyEnd(xvec,ierr)
@ -99,34 +103,6 @@ print *,'solving system'
call VecRestoreArrayF90(xvec,xptr,ierr)
end if
! allocate and create number of nonzeros vector
if (row_end /= n) then
if (.not. allocated(nnzv)) allocate(nnzv(row_start:row_end-1))
nnzv = loss%nnz + 1
else
if (.not. allocated(nnzv)) allocate(nnzv(row_start:row_end))
nnzv = loss%nnz + 1
nnzv(n) = n + 1
end if
! get ownership rows again
call VecGetOwnershipRange(xvec,row_start,row_end,ierr)
! set local size
n = row_end - row_start
! take minimum of size (incase the local size is less than the number of nnz)
nnzv = min(n,nnzv)
! create the preconditioner matrix
call MatCreateMPIAIJ(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,PETSC_NULL_INTEGER, &
& nnzv,PETSC_NULL_INTEGER,nnzv,jac_prec,ierr)
call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr)
call MatSetOption(jac_prec,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr)
call MatGetOwnershipRange(jac_prec,row_start,row_end,ierr)
call MPI_Barrier(MPI_COMM_WORLD,ierr)
stop
end subroutine init_data
!===============================================================================
@ -161,7 +137,7 @@ stop
call MatCreateSNESMF(snes,jac,ierr)
! set matrix free finite difference
call SNESSetJacobian(snes,jac,jac_prec,compute_jacobian,PETSC_NULL,ierr)
call SNESSetJacobian(snes,jac,jac_prec,build_jacobian_matrix,ctx,ierr)
! set lags
call SNESSetLagJacobian(snes,-2,ierr)
@ -209,15 +185,20 @@ stop
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
integer :: rank ! rank of processor
integer :: n_procs ! number of processors
real(8) :: lambda ! eigenvalue
real(8) :: reslamb ! residual for lambda
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: rptr(:) ! pointer to residual vector
print *,'In residual'
! get number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD,n_procs,ierr)
! get problem size
n = loss%n
call VecGetSize(x,n,ierr)
n = n - 1 ! subtract off last row
! get the local dimensions for each process
call VecGetOwnershipRange(x,row_start,row_end,ierr)
@ -234,31 +215,29 @@ stop
call VecPlaceArray(phi,xptr,ierr)
call VecPlaceArray(rphi,rptr,ierr)
! extract eigenvalue and broadcast
! extract eigenvalue and broadcast (going to want to make this more general in future)
if (row_end == n+1) then
lambda = xptr(n+1)
call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
call MPI_BCAST(lambda,1,MPI_REAL8,rank,MPI_COMM_WORLD,ierr)
lambda = xptr(size(xptr))
end if
print *,lambda
stop
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
! create operators
call build_loss_matrix(loss)
call build_prod_matrix(prod)
call build_loss_matrix(ctx%loss)
call build_prod_matrix(ctx%prod)
! create new petsc vectors to perform math
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,phiM,ierr)
! calculate flux part of residual vector
call MatMult(loss%M,phi,phiM,ierr)
call MatMult(prod%F,phi,rphi,ierr)
call MatMult(ctx%loss%M,phi,phiM,ierr)
call MatMult(ctx%prod%F,phi,rphi,ierr)
call VecAYPX(rphi,-1.0_8*lambda,phiM,ierr)
! set eigenvalue part of residual vector
call VecDot(phi,phi,reslamb,ierr)
! map to ptr
if (row_end == n+1) rptr(n+1) = 0.5_8 - 0.5_8*reslamb
if (row_end == n+1) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
! reset arrays that are not used
call VecResetArray(phi,ierr)
@ -275,119 +254,6 @@ stop
end subroutine compute_nonlinear_residual
!===============================================================================
! COMPUTE_JACOBIAN
!===============================================================================
subroutine compute_jacobian(snes,x,jac,jac_prec,flag,user,ierr)
! formal variables
SNES :: snes ! the snes context
Vec :: x ! the solution vector
Mat :: jac ! the jacobian matrix
Mat :: jac_prec ! the jacobian preconditioner
MatStructure :: flag ! not used
integer :: user(*) ! not used
integer :: ierr ! petsc error flag
! local varibles
Vec :: phi ! flux vector
Vec :: source ! source vector
integer :: n ! problem size
integer :: k ! implied do loop counter
integer :: ncols ! number of nonzeros in cols
integer :: irow ! row counter
integer, allocatable :: nnzv(:) ! vector of number of nonzeros for jac
integer, allocatable :: cols(:) ! vector of column numbers
real(8) :: lambda ! eigenvalue
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: sptr(:) ! pointer to source vector
real(8), allocatable :: vals(:) ! vector of row values
! create operators
call build_loss_matrix(loss)
call build_prod_matrix(prod)
! get problem size
n = loss%n
! allocate cols and rows and initialize to zero
if (.not. allocated(cols)) allocate(cols(loss%nnz))
if (.not. allocated(vals)) allocate(vals(loss%nnz))
cols = 0
vals = 0.0_8
! get pointers to residual vector
call VecGetArrayF90(x,xptr,ierr)
! create petsc vector for flux
call VecCreate(PETSC_COMM_SELF,phi,ierr)
call VecSetSizes(phi,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(phi,ierr)
! extract flux and eigenvalue
call VecPlaceArray(phi,xptr,ierr)
lambda = xptr(n+1)
! compute math (M-lambda*F) M is overwritten here
call MatAXPY(loss%M,-1.0_8*lambda,prod%F,SUBSET_NONZERO_PATTERN,ierr)
! create tmp petsc vector for source
call VecCreate(PETSC_COMM_SELF,source,ierr)
call VecSetSizes(source,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(source,ierr)
! perform math (-F*phi --> source)
call MatMult(prod%F,phi,source,ierr)
call VecScale(source,-1.0_8,ierr)
! get pointer to source
call VecGetArrayF90(source,sptr,ierr)
! begin loop to insert things into matrix
do irow = 0,n-1
! get row of matrix
call MatGetRow(loss%M,irow,ncols,cols,vals,ierr)
! set that row to Jacobian matrix
call MatSetValues(jac_prec,1,irow,ncols,cols(1:ncols),vals,INSERT_VALUES,ierr)
! restore the row
call MatRestoreRow(loss%M,irow,ncols,cols,vals,ierr)
! insert source value
call MatSetValue(jac_prec,irow,n,sptr(irow+1),INSERT_VALUES,ierr)
end do
! set values in last row of matrix
call MatSetValues(jac_prec,1,n,n,(/(k,k=0,n-1)/),xptr(1:n),INSERT_VALUES,ierr)
call MatSetValue(jac_prec,n,n,1.0_8,INSERT_VALUES,ierr)
! assemble matrix
call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
! reset all vectors
call VecResetArray(phi,ierr)
! restore all vectors
call VecRestoreArrayF90(x,xptr,ierr)
call VecRestoreArrayF90(source,sptr,ierr)
! destroy all temporary objects
call VecDestroy(phi,ierr)
call VecDestroy(source,ierr)
! deallocate all temporary space
if (allocated(cols)) deallocate(cols)
if (allocated(vals)) deallocate(vals)
end subroutine compute_jacobian
!===============================================================================
! EXTRACT_RESULTS
!===============================================================================
@ -400,7 +266,7 @@ stop
PetscScalar, pointer :: xptr(:) ! pointer to eigenvector info
! get problem size
n = loss%n
n = ctx%loss%n
! also allocate in cmfd object
if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n))
@ -422,8 +288,9 @@ stop
subroutine finalize()
! finalize data objects
call destroy_M_operator(loss)
call destroy_F_operator(prod)
call destroy_M_operator(ctx%loss)
call destroy_F_operator(ctx%prod)
call destroy_J_operator(jac_prec)
call VecDestroy(xvec,ierr)
call VecDestroy(resvec,ierr)