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 listsmodule, only: basesolutionlist
  use linearsolverbasemodule
  use matrixbasemodule
  use vectorbasemodule
  use petscmatrixmodule
  use petscvectormodule
  use petscconvergencemodule
  use petscimspreconditionermodule
  use convergencesummarymodule
  use imslinearsettingsmodule
  use simvariablesmodule, only: iout, simulation_mode, nr_procs, proc_id
  use simmodule, only: store_error, store_warning
  use devfeaturemodule, only: dev_feature
 
  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
    logical(LGP) :: use_ims_cnvgopt !< when true, use IMS convergence check in PETSc solve
    ksptype :: ksp_type !< the KSP solver type (CG, BCGS, ...)
    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
    character(len=LENSOLUTIONNAME + 1) :: option_prefix !< prefix for keys in petscrc database in case there are multiple solutions
  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_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)
    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
 
  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_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)
 
  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
    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
 
  end subroutine get_options_mf6
 
  !> @brief Create the PETSc KSP object
  !<
  subroutine create_ksp(this)
    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)
 
  end subroutine create_ksp
 
  !> @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
 
    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)
 
  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)
 
    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)
 
  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
    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
 
  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