PetscSolver.F90

Go to the documentation of this file.

module petscsolvermodule
#include <petsc/finclude/petscksp.h>
  use petscksp
  use kindmodule, only: i4b, dp, lgp
  use constantsmodule, only: linelength, lensolutionname, dzero
  use linearsolverbasemodule
  use matrixbasemodule
  use vectorbasemodule
  use petscmatrixmodule
  use petscvectormodule
  use petscconvergencemodule
  use petscimspreconditionermodule
  use convergencesummarymodule
  use imslinearsettingsmodule
  use simvariablesmodule, only: iout, simulation_mode
  use simmodule, only: store_error, store_warning
 
  implicit none
  private
 
  public :: create_petsc_solver
 
  !< Linear solver using PETSc KSP
  type, public, extends(linearsolverbasetype) :: petscsolvertype
    ksp :: ksp_petsc !< the KSP linear solver object
    class(petscmatrixtype), pointer :: matrix => null() !< the system matrix in PETSc compatible format
    mat, pointer :: mat_petsc => null() !< the PETSc internal matrix format
 
    type(imslinearsettingstype), pointer :: linear_settings => null() !< pointer to linear settings from IMS
    logical(LGP) :: use_ims_pc !< when true, use custom IMS-style preconditioning
    ksptype :: ksp_type !< the KSP solver type (CG, BCGS, ...)
    pctype :: pc_type !< the (global) preconditioner type, should be parallel
    pctype :: sub_pc_type !< the block preconditioner type (for the subdomain)
    class(petsccnvgctxtype), pointer :: petsc_ctx => null() !< context for the PETSc custom convergence check
    type(pcshellctxtype), pointer :: pc_context => null() !< context for the custom (IMS) precondioner
    type(convergencesummarytype), pointer :: convergence_summary => null() !< data structure wrapping the convergence data
 
  contains
    procedure :: initialize => petsc_initialize
    procedure :: solve => petsc_solve
    procedure :: print_summary => petsc_print_summary
    procedure :: destroy => petsc_destroy
    procedure :: create_matrix => petsc_create_matrix
 
    ! private
    procedure, private :: petsc_check_settings
    procedure, private :: get_options_mf6
    procedure, private :: create_ksp
    procedure, private :: create_convergence_check
    procedure, private :: set_petsc_pc
    procedure, private :: set_ims_pc
    procedure, private :: print_vec
    procedure, private :: print_petsc_version
  end type petscsolvertype
 
contains
 
  !> @brief Create a PETSc solver object
  !<
  function create_petsc_solver(sln_name) result(solver)
    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
 
  end function create_petsc_solver
 
  !> @brief Initialize PETSc KSP solver with
  !<  options from the petsc database file
  subroutine petsc_initialize(this, matrix, linear_settings, convergence_summary)
    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)
 
  end subroutine petsc_initialize
 
  subroutine petsc_check_settings(this, linear_settings)
    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
 
  end subroutine petsc_check_settings
 
  !> @brief Print PETSc version string from shared lib
  !<
  subroutine print_petsc_version(this)
    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
 
  end subroutine print_petsc_version
 
  !> @brief Get the MODFLOW specific options from the PETSc database
  !<
  subroutine get_options_mf6(this)
    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)
 
  end subroutine get_options_mf6
 
  !> @brief Create the PETSc KSP object
  !<
  subroutine create_ksp(this)
    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)
 
  end subroutine create_ksp
 
  !> @brief Set up a standard PETSc preconditioner from
  !< the configured settings
  subroutine set_petsc_pc(this)
    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)
 
  end subroutine set_petsc_pc
 
  !> @brief Set up a custom preconditioner following the ones
  !< we have in IMS, i.e. Modified ILU(T)
  subroutine set_ims_pc(this)
    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)
 
  end subroutine set_ims_pc
 
  !> @brief Create and assign a custom convergence
  !< check for this solver
  subroutine create_convergence_check(this, convergence_summary)
    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)
 
  end subroutine create_convergence_check
 
  subroutine petsc_solve(this, kiter, rhs, x, cnvg_summary)
    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
 
  end subroutine petsc_solve
 
  subroutine petsc_print_summary(this)
    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))
 
  end subroutine petsc_print_summary
 
  subroutine petsc_destroy(this)
    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)
 
  end subroutine petsc_destroy
 
  function petsc_create_matrix(this) result(matrix)
    class(petscsolvertype) :: this
    class(matrixbasetype), pointer :: matrix
    ! local
    class(petscmatrixtype), pointer :: petsc_matrix
 
    allocate (petsc_matrix)
    matrix => petsc_matrix
 
  end function petsc_create_matrix
 
  subroutine print_vec(this, vec, vec_name, kiter)
    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)
 
  end subroutine print_vec
 
end module petscsolvermodule