diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 166e71074a..30fbda4c3e 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/Makefile b/src/Makefile index 0438f1bedb..afc6b8ab5f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index f79027e655..1a9f2b8cc7 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -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 \ diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 008d299708..f1905d34ce 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -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() diff --git a/src/cmfd_jacobian_operator.F90 b/src/cmfd_jacobian_operator.F90 new file mode 100644 index 0000000000..745b59fc88 --- /dev/null +++ b/src/cmfd_jacobian_operator.F90 @@ -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 + + 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 diff --git a/src/cmfd_slepc_solver.F90 b/src/cmfd_slepc_solver.F90 index 5e903e000a..9fcffedf93 100644 --- a/src/cmfd_slepc_solver.F90 +++ b/src/cmfd_slepc_solver.F90 @@ -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() diff --git a/src/cmfd_snes_solver.F90 b/src/cmfd_snes_solver.F90 index 8d817bd72b..9cd8a7d232 100644 --- a/src/cmfd_snes_solver.F90 +++ b/src/cmfd_snes_solver.F90 @@ -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 - 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)