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_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 291 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)
 
    if (.not. this%use_ims_cnvgopt) then
      ! use PETSc residual L2 norm for convergence
      call dev_feature('Using PETSc convergence is under development, install &
      &the nightly build or compile from source with IDEVELOPMODE = 1.')
      call kspsetconvergencetest(this%ksp_petsc, petsc_cnvg_check_internal, &
                                 this%petsc_ctx, petsc_null_function, ierr)
    else
      ! IMS convergence check
      call kspsetconvergencetest(this%ksp_petsc, petsc_cnvg_check, &
                                 this%petsc_ctx, petsc_null_function, ierr)
    end if
    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 210 of file PetscSolver.F90.

    class(PetscSolverType) :: this !< This solver instance
    ! local
    petscerrorcode :: ierr
    pc :: pc
 
    call kspcreate(petsc_comm_world, this%ksp_petsc, ierr)
    chkerrq(ierr)
 
    ! set prefix for options database
    call kspsetoptionsprefix(this%ksp_petsc, trim(this%option_prefix), 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 dev_feature('PETSc preconditioning is under development, install the &
            &nightly build or compile from source with IDEVELOPMODE = 1.')
      ! The PC options will be set from the .petscrc
      ! file in the call to KSPSetFromOptions below
      call kspgetpc(this%ksp_petsc, pc, ierr)
      chkerrq(ierr)
      call pcsetfromoptions(pc, ierr)
      chkerrq(ierr)
    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)
 
    call kspsetup(this%ksp_petsc, ierr)
    chkerrq(ierr)
 

Here is the call graph for this function:

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)
    petsc_solver%option_prefix = ""
    if (basesolutionlist%Count() > 1) then
      petsc_solver%option_prefix = trim(sln_name)//"_"
    end if
 
    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 187 of file PetscSolver.F90.

    class(PetscSolverType) :: this
    ! local
    petscerrorcode :: ierr
    logical(LGP) :: found
    logical(LGP) :: use_petsc_pc, use_petsc_cnvg
 
    use_petsc_pc = .false.
    call petscoptionsgetbool(petsc_null_options, trim(this%option_prefix), &
                             '-use_petsc_pc', use_petsc_pc, found, ierr)
    chkerrq(ierr)
    this%use_ims_pc = .not. use_petsc_pc
 
    use_petsc_cnvg = .false.
    call petscoptionsgetbool(petsc_null_options, trim(this%option_prefix), &
                             '-use_petsc_cnvg', use_petsc_cnvg, found, ierr)
    chkerrq(ierr)
    this%use_ims_cnvgopt = .not. use_petsc_cnvg
 

petsc_check_settings()

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

private

Definition at line 128 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 468 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 451 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 87 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_cnvgopt = .true. ! use IMS convergence check, override with .petscrc
    this%use_ims_pc = .true. ! use IMS preconditioning, 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
 
    ! 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 384 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
    character(len=128) :: ksp_logfile
    integer :: ierr
    pc :: pc
    petscviewer :: ksp_viewer
 
    if (this%use_ims_pc) then
      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)
      subpc_str = this%pc_context%ims_pc_type
 
      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, '(/,1x,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))
      if (this%use_ims_cnvgopt) then
        write (iout, '(1x,a,i0)') &
          "Residual convergence option:  ", this%linear_settings%icnvgopt
      else
        write (iout, '(1x,a)') &
          "Residual convergence option:  PETSc L2 norm"
      end if
      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))
    else
      ksp_logfile = trim(this%option_prefix)//"ksp_logview.txt"
      write (iout, '(/,1x,a)') "PETSc linear solver settings from .petscrc: "
      write (iout, '(1x,a)') repeat('-', 66)
      write (iout, '(1x,2a)') "see ", trim(ksp_logfile)
 
      ! collective write
      call petscviewerasciiopen(petsc_comm_world, ksp_logfile, ksp_viewer, ierr); 
      chkerrq(ierr)
      call petscviewerpushformat(ksp_viewer, petsc_viewer_ascii_info_detail, ierr)
      chkerrq(ierr)
      call kspview(this%ksp_petsc, ksp_viewer, ierr)
      chkerrq(ierr)
      call petscviewerdestroy(ksp_viewer, ierr)
      chkerrq(ierr)
    end if
 

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 316 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) then
      if (cnvg_reason == ksp_diverged_breakdown .or. &
          cnvg_reason == ksp_diverged_its) then
        ! out of iterations, or rho/eta became zero,
        ! move to next Picard iteration and try again
        this%is_converged = 0
      else
        ! this is an solver error, terminate
        write (errmsg, '(1x,3a,i0)') "PETSc convergence failure in ", &
          trim(this%name), ": ", cnvg_reason
        call store_error(errmsg, terminate=.true.)
      end if
    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 163 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 479 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 260 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, pcbjacobi, 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, pcshell, 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: