petscsolvermodule Module Reference

Data Types
type	petscsolvertype

Functions/Subroutines
class(linearsolverbasetype) function, pointer, public	create_petsc_solver (sln_name)
	Create a PETSc solver object. More...
subroutine	petsc_initialize (this, matrix, linear_settings, convergence_summary)
	Initialize PETSc KSP solver with. More...
subroutine	petsc_check_settings (this, linear_settings)
subroutine	print_petsc_version (this)
	Print PETSc version string from shared lib. More...
subroutine	get_options_mf6 (this)
	Get the MODFLOW specific options from the PETSc database. More...
subroutine	create_ksp (this)
	Create the PETSc KSP object. More...
subroutine	set_petsc_pc (this)
	Set up a standard PETSc preconditioner from. More...
subroutine	set_ims_pc (this)
	Set up a custom preconditioner following the ones. More...
subroutine	create_convergence_check (this, convergence_summary)
	Create and assign a custom convergence. More...
subroutine	petsc_solve (this, kiter, rhs, x, cnvg_summary)
subroutine	petsc_print_summary (this)
subroutine	petsc_destroy (this)
class(matrixbasetype) function, pointer	petsc_create_matrix (this)
subroutine	print_vec (this, vec, vec_name, kiter)

Function/Subroutine Documentation

create_convergence_check()

subroutine petscsolvermodule::create_convergence_check ( class(petscsolvertype)  this, type(convergencesummarytype), pointer  convergence_summary  )

private

Parameters

this	This solver instance

Definition at line 296 of file PetscSolver.F90.

    use imslinearbasemodule, only: ims_base_epfact
    class(PetscSolverType) :: this !< This solver instance
    type(ConvergenceSummaryType), pointer :: convergence_summary
    ! local
    petscerrorcode :: ierr
 
    call this%petsc_ctx%create(this%mat_petsc, this%linear_settings, &
                               convergence_summary)
 
    call kspsetconvergencetest(this%ksp_petsc, petsc_check_convergence, &
                               this%petsc_ctx, petsc_null_function, ierr)
    chkerrq(ierr)
 

Here is the call graph for this function:

create_ksp()

subroutine petscsolvermodule::create_ksp ( class(petscsolvertype)  this )

private

Parameters

this	This solver instance

Definition at line 199 of file PetscSolver.F90.

    class(PetscSolverType) :: this !< This solver instance
    ! local
    petscerrorcode :: ierr
 
    call kspcreate(petsc_comm_world, this%ksp_petsc, ierr)
    chkerrq(ierr)
 
    call kspsetoperators(this%ksp_petsc, this%mat_petsc, this%mat_petsc, ierr)
    chkerrq(ierr)
 
    call kspsetinitialguessnonzero(this%ksp_petsc, .true., ierr)
    chkerrq(ierr)
 
    call kspsettype(this%ksp_petsc, this%ksp_type, ierr)
    chkerrq(ierr)
 
    if (this%use_ims_pc) then
      call this%set_ims_pc()
    else
      call this%set_petsc_pc()
    end if
 
    call kspseterrorifnotconverged(this%ksp_petsc, .false., ierr)
    chkerrq(ierr)
 
    ! finally override these options from the
    ! optional .petscrc file
    call kspsetfromoptions(this%ksp_petsc, ierr)
    chkerrq(ierr)
 

create_petsc_solver()

class(linearsolverbasetype) function, pointer, public petscsolvermodule::create_petsc_solver

(

character(len=lensolutionname)

sln_name

)

Returns

Uninitialized instance of the PETSc solver

Parameters

sln_name

the solution name

Definition at line 60 of file PetscSolver.F90.

    class(LinearSolverBaseType), pointer :: solver !< Uninitialized instance of the PETSc solver
    character(len=LENSOLUTIONNAME) :: sln_name !< the solution name
    ! local
    class(PetscSolverType), pointer :: petsc_solver
    character(len=LINELENGTH) :: errmsg
 
    allocate (petsc_solver)
    allocate (petsc_solver%petsc_ctx)
 
    solver => petsc_solver
    solver%name = sln_name
 
    if (simulation_mode /= 'PARALLEL') then
      write (errmsg, '(a,a)') 'PETSc solver not supported for run mode: ', &
        trim(simulation_mode)
      call store_error(errmsg, terminate=.true.)
    end if
 

Here is the call graph for this function:

Here is the caller graph for this function:

get_options_mf6()

subroutine petscsolvermodule::get_options_mf6 ( class(petscsolvertype)  this )

private

Definition at line 185 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    ! local
    petscerrorcode :: ierr
    logical(LGP) :: found
 
    call petscoptionsgetbool(petsc_null_options, petsc_null_character, &
                             '-ims_pc', this%use_ims_pc, found, ierr)
    chkerrq(ierr)
 

petsc_check_settings()

subroutine petscsolvermodule::petsc_check_settings ( class(petscsolvertype)  this, type(imslinearsettingstype), pointer  linear_settings  )

private

Definition at line 126 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    type(ImsLinearSettingsType), pointer :: linear_settings
    ! local
    character(len=LINELENGTH) :: warnmsg, errmsg
 
    ! errors
    if (linear_settings%ilinmeth /= cg_method .and. &
        linear_settings%ilinmeth /= bcgs_method) then
      write (errmsg, '(a,a)') 'PETSc: unknown linear solver method in ', &
        this%name
      call store_error(errmsg, terminate=.true.)
    end if
 
    ! warnings
    if (linear_settings%iord > 0) then
      linear_settings%iord = 0
      write (warnmsg, '(a)') 'PETSc: IMS reordering not supported'
      call store_warning(warnmsg)
    end if
    if (linear_settings%iscl > 0) then
      linear_settings%iscl = 0
      write (warnmsg, '(a)') 'PETSc: IMS matrix scaling not supported'
      call store_warning(warnmsg)
    end if
    if (linear_settings%north > 0) then
      linear_settings%north = 0
      write (warnmsg, '(a)') 'PETSc: IMS orthogonalization not supported'
      call store_warning(warnmsg)
    end if
 

Here is the call graph for this function:

petsc_create_matrix()

class(matrixbasetype) function, pointer petscsolvermodule::petsc_create_matrix ( class(petscsolvertype)  this )

private

Definition at line 436 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    class(MatrixBaseType), pointer :: matrix
    ! local
    class(PetscMatrixType), pointer :: petsc_matrix
 
    allocate (petsc_matrix)
    matrix => petsc_matrix
 

petsc_destroy()

subroutine petscsolvermodule::petsc_destroy ( class(petscsolvertype)  this )

private

Definition at line 419 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    ! local
    petscerrorcode :: ierr
 
    call kspdestroy(this%ksp_petsc, ierr)
    chkerrq(ierr)
 
    ! delete context
    call this%petsc_ctx%destroy()
    deallocate (this%petsc_ctx)
 
    call this%pc_context%destroy()
    deallocate (this%pc_context)
 

petsc_initialize()

subroutine petscsolvermodule::petsc_initialize ( class(petscsolvertype)  this, class(matrixbasetype), pointer  matrix, type(imslinearsettingstype), pointer  linear_settings, type(convergencesummarytype), pointer  convergence_summary  )

private

Parameters

this	This solver instance
matrix	The solution matrix as KSP operator
linear_settings	the settings for the linear solver from the .ims file
convergence_summary	a convergence record for diagnostics

Definition at line 83 of file PetscSolver.F90.

    class(PetscSolverType) :: this !< This solver instance
    class(MatrixBaseType), pointer :: matrix !< The solution matrix as KSP operator
    type(ImsLinearSettingsType), pointer :: linear_settings !< the settings for the linear solver from the .ims file
    type(ConvergenceSummaryType), pointer :: convergence_summary !< a convergence record for diagnostics
 
    this%linear_settings => linear_settings
 
    call this%print_petsc_version()
 
    this%mat_petsc => null()
    select type (pm => matrix)
    class is (petscmatrixtype)
      this%matrix => pm
      this%mat_petsc => pm%mat
    end select
 
    call this%petsc_check_settings(linear_settings)
 
    this%use_ims_pc = .true. ! default is true, override with .petscrc
    allocate (this%pc_context)
    call this%pc_context%create(this%matrix, linear_settings)
 
    if (linear_settings%ilinmeth == cg_method) then
      this%ksp_type = kspcg
    else
      this%ksp_type = kspbcgs
    end if
 
    this%pc_type = pcbjacobi
    this%sub_pc_type = pcilu
 
    ! get MODFLOW options from PETSc database file
    call this%get_options_mf6()
 
    ! create the solver object
    call this%create_ksp()
 
    ! Create custom convergence check
    call this%create_convergence_check(convergence_summary)
 

petsc_print_summary()

subroutine petscsolvermodule::petsc_print_summary ( class(petscsolvertype)  this )

private

Definition at line 372 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    ! local
    character(len=128) :: ksp_str, pc_str, subpc_str, &
                          dvclose_str, rclose_str, relax_str, dtol_str
    integer :: ierr
    pc :: pc
 
    call kspgettype(this%ksp_petsc, ksp_str, ierr)
    chkerrq(ierr)
    call kspgetpc(this%ksp_petsc, pc, ierr)
    chkerrq(ierr)
    call pcgettype(pc, pc_str, ierr)
    chkerrq(ierr)
    if (this%use_ims_pc) then
      subpc_str = this%pc_context%ims_pc_type
    else
      subpc_str = this%sub_pc_type
    end if
 
    write (dvclose_str, '(e15.5)') this%linear_settings%dvclose
    write (rclose_str, '(e15.5)') this%linear_settings%rclose
    write (relax_str, '(e15.5)') this%linear_settings%relax
    write (dtol_str, '(e15.5)') this%linear_settings%droptol
 
    write (iout, '(/,7x,a)') "PETSc linear solver settings: "
    write (iout, '(1x,a)') repeat('-', 66)
    write (iout, '(1x,a,a)') "Linear acceleration method:   ", trim(ksp_str)
    write (iout, '(1x,a,a)') "Preconditioner type:          ", trim(pc_str)
    write (iout, '(1x,a,a)') "Sub-preconditioner type:      ", trim(subpc_str)
    write (iout, '(1x,a,i0)') "Maximum nr. of iterations:    ", &
      this%linear_settings%iter1
    write (iout, '(1x,a,a)') &
      "Dep. var. closure criterion:  ", trim(adjustl(dvclose_str))
    write (iout, '(1x,a,a)') &
      "Residual closure criterion:   ", trim(adjustl(rclose_str))
    write (iout, '(1x,a,i0)') &
      "Residual convergence option:  ", this%linear_settings%icnvgopt
    write (iout, '(1x,a,a)') &
      "Relaxation factor MILU(T):    ", trim(adjustl(relax_str))
    write (iout, '(1x,a,i0)') &
      "Fill level in factorization:  ", this%linear_settings%level
    write (iout, '(1x,a,a,/)') &
      "Drop tolerance level fill:    ", trim(adjustl(dtol_str))
 

petsc_solve()

subroutine petscsolvermodule::petsc_solve	(	class(petscsolvertype)	this,
		integer(i4b)	kiter,
		class(vectorbasetype), pointer	rhs,
		class(vectorbasetype), pointer	x,
		type(convergencesummarytype)	cnvg_summary
	)

Definition at line 312 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    integer(I4B) :: kiter
    class(VectorBaseType), pointer :: rhs
    class(VectorBaseType), pointer :: x
    type(ConvergenceSummaryType) :: cnvg_summary
    ! local
    petscerrorcode :: ierr
    class(PetscVectorType), pointer :: rhs_petsc, x_petsc
    kspconvergedreason :: cnvg_reason
    integer :: it_number
    character(len=LINELENGTH) :: errmsg
 
    rhs_petsc => null()
    select type (rhs)
    class is (petscvectortype)
      rhs_petsc => rhs
    end select
 
    x_petsc => null()
    select type (x)
    class is (petscvectortype)
      x_petsc => x
    end select
 
    this%iteration_number = 0
    this%is_converged = 0
    if (kiter == 1) then
      this%petsc_ctx%cnvg_summary%iter_cnt = 0
    end if
 
    ! update matrix coefficients
    call this%matrix%update()
    call kspsolve(this%ksp_petsc, rhs_petsc%vec_impl, x_petsc%vec_impl, ierr)
    chkerrq(ierr)
 
    call kspgetiterationnumber(this%ksp_petsc, it_number, ierr)
    chkerrq(ierr)
    this%iteration_number = it_number
 
    call kspgetconvergedreason(this%ksp_petsc, cnvg_reason, ierr)
    chkerrq(ierr)
    if (cnvg_reason > 0) then
      if (this%petsc_ctx%icnvg_ims == -1) then
        ! move to next Picard iteration (e.g. with 'STRICT' option)
        this%is_converged = 0
      else
        ! linear convergence reached
        this%is_converged = 1
      end if
    end if
 
    if (cnvg_reason < 0 .and. cnvg_reason /= ksp_diverged_its) then
      write (errmsg, '(1x,3a,i0)') "PETSc convergence failure in ", &
        trim(this%name), ": ", cnvg_reason
      call store_error(errmsg, terminate=.true.)
    end if
 

Here is the call graph for this function:

print_petsc_version()

subroutine petscsolvermodule::print_petsc_version ( class(petscsolvertype)  this )

private

Definition at line 161 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    ! local
    petscerrorcode :: ierr
    petscint :: major, minor, subminor, release
    character(len=128) :: petsc_version, release_str
 
    call petscgetversionnumber(major, minor, subminor, release, ierr)
    chkerrq(ierr)
 
    if (release == 1) then
      release_str = "(release)"
    else
      release_str = "(unofficial)"
    end if
    write (petsc_version, '(i0,a,i0,a,i0,a,a)') &
      major, ".", minor, ".", subminor, " ", trim(release_str)
    write (iout, '(/,1x,4a,/)') "PETSc Linear Solver will be used for ", &
      trim(this%name), ": version ", petsc_version
 

print_vec()

subroutine petscsolvermodule::print_vec ( class(petscsolvertype)  this, class(petscvectortype)  vec, character(len=*)  vec_name, integer(i4b)  kiter  )

private

Definition at line 447 of file PetscSolver.F90.

    use tdismodule, only: nper, kstp
    class(PetscSolverType) :: this
    class(PetscVectorType) :: vec
    character(len=*) :: vec_name
    integer(I4B) :: kiter
    ! local
    petscviewer :: viewer
    character(len=24) :: filename
    petscerrorcode :: ierr
 
    write (filename, '(2a,i0,a,i0,a,i0,a)') vec_name, '_', nper, &
      '_', kstp, '_', kiter, '.txt'
    call petscviewerasciiopen(petsc_comm_world, filename, viewer, ierr)
    chkerrq(ierr)
    call vecview(vec%vec_impl, viewer, ierr)
    chkerrq(ierr)
    call petscviewerdestroy(viewer, ierr)
    chkerrq(ierr)
 

set_ims_pc()

subroutine petscsolvermodule::set_ims_pc ( class(petscsolvertype)  this )

private

Parameters

this	This solver instance

Definition at line 263 of file PetscSolver.F90.

    class(PetscSolverType) :: this !< This solver instance
    ! local
    pc :: pc, sub_pc
    ksp, dimension(1) :: sub_ksp
    petscint :: n_local, n_first
    petscerrorcode :: ierr
 
    this%sub_pc_type = pcshell
 
    call kspgetpc(this%ksp_petsc, pc, ierr)
    chkerrq(ierr)
    call pcsettype(pc, this%pc_type, ierr)
    chkerrq(ierr)
    call pcsetup(pc, ierr)
    chkerrq(ierr)
    call pcbjacobigetsubksp(pc, n_local, n_first, sub_ksp, ierr)
    chkerrq(ierr)
    call kspgetpc(sub_ksp(1), sub_pc, ierr)
    chkerrq(ierr)
    call pcsettype(sub_pc, this%sub_pc_type, ierr)
    chkerrq(ierr)
    call pcshellsetapply(sub_pc, pcshell_apply, ierr)
    chkerrq(ierr)
    call pcshellsetsetup(sub_pc, pcshell_setup, ierr)
    chkerrq(ierr)
    call pcshellsetcontext(sub_pc, this%pc_context, ierr)
    chkerrq(ierr)
 

Here is the call graph for this function:

set_petsc_pc()

subroutine petscsolvermodule::set_petsc_pc ( class(petscsolvertype)  this )

private

Parameters

this	This solver instance

Definition at line 234 of file PetscSolver.F90.

    class(PetscSolverType) :: this !< This solver instance
    ! local
    pc :: pc, sub_pc
    ksp, dimension(1) :: sub_ksp
    petscint :: n_local, n_first
    petscerrorcode :: ierr
 
    call kspgetpc(this%ksp_petsc, pc, ierr)
    chkerrq(ierr)
    call pcsettype(pc, this%pc_type, ierr)
    chkerrq(ierr)
    call pcsetup(pc, ierr)
    chkerrq(ierr)
    call pcbjacobigetsubksp(pc, n_local, n_first, sub_ksp, ierr)
    chkerrq(ierr)
    call kspgetpc(sub_ksp(1), sub_pc, ierr)
    chkerrq(ierr)
    call pcsettype(sub_pc, this%sub_pc_type, ierr)
    chkerrq(ierr)
    call pcsetfromoptions(sub_pc, ierr)
    chkerrq(ierr)
    call pcfactorsetlevels(sub_pc, this%linear_settings%level, ierr)
    chkerrq(ierr)