Actual source code: ex45f.F90

      1: #include <petsc/finclude/petscdmda.h>
      2: #include <petsc/finclude/petscksp.h>
      3: module ex45fmodule
      4:   use petscmat
      5:   use petscksp
      6:   use petscdmda
      7:   use petscvec
      8:   implicit none

     10: contains

     12:   subroutine ComputeInitialGuess(ksp, b, ctx, ierr)
     13:     PetscErrorCode, intent(out) :: ierr
     14:     KSP ksp
     15:     PetscInt ctx(*)
     16:     Vec b
     17:     PetscScalar, parameter :: h = 0.0

     19:     PetscCall(VecSet(b, h, ierr))
     20:   end subroutine

     22:   subroutine ComputeRHS(ksp, b, unused, ierr)
     23:     PetscErrorCode, intent(out) :: ierr
     24:     KSP ksp
     25:     Vec b
     26:     integer unused(*)
     27:     PetscScalar h, Hx, Hy
     28:     PetscInt mx, my
     29:     DM dm

     31:     PetscCall(KSPGetDM(ksp, dm, ierr))
     32:     PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

     34:     Hx = 1.0/real(mx - 1)
     35:     Hy = 1.0/real(my - 1)
     36:     h = Hx*Hy
     37:     PetscCall(VecSet(b, h, ierr))
     38:   end subroutine

     40:   subroutine ComputeMatrix(ksp, A, B, unused, ierr)
     41:     PetscErrorCode, intent(out) :: ierr
     42:     KSP ksp
     43:     Mat A, B
     44:     integer unused(*)
     45:     DM dm

     47:     PetscInt i, j, mx, my, xm, ym, xs, ys
     48:     PetscScalar v(5), Hx, Hy
     49:     PetscScalar HxdHy, HydHx
     50:     MatStencil row(1), col(5)

     52:     PetscCall(KSPGetDM(ksp, dm, ierr))
     53:     PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

     55:     Hx = 1.0/real(mx - 1)
     56:     Hy = 1.0/real(my - 1)
     57:     HxdHy = Hx/Hy
     58:     HydHx = Hy/Hx
     59:     PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
     60:     do j = ys, ys + ym - 1
     61:       do i = xs, xs + xm - 1
     62:         row(1)%i = i
     63:         row(1)%j = j
     64:         if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
     65:           v(1) = 2.0*(HxdHy + HydHx)
     66:           PetscCall(MatSetValuesStencil(B, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, row, v, INSERT_VALUES, ierr))
     67:         else
     68:           v(1) = -HxdHy
     69:           col(1)%i = i
     70:           col(1)%j = j - 1
     71:           v(2) = -HydHx
     72:           col(2)%i = i - 1
     73:           col(2)%j = j
     74:           v(3) = 2.0*(HxdHy + HydHx)
     75:           col(3)%i = i
     76:           col(3)%j = j
     77:           v(4) = -HydHx
     78:           col(4)%i = i + 1
     79:           col(4)%j = j
     80:           v(5) = -HxdHy
     81:           col(5)%i = i
     82:           col(5)%j = j + 1
     83:           PetscCall(MatSetValuesStencil(B, 1_PETSC_INT_KIND, row, 5_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
     84:         end if
     85:       end do
     86:     end do
     87:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
     88:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
     89:     if (A /= B) then
     90:       PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
     91:       PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
     92:     end if
     93:   end subroutine
     94: end module ex45fmodule

     96: program main
     97:   use petscdmda
     98:   use petscksp
     99:   use ex45fmodule
    100:   implicit none

    102:   PetscInt is, js, iw, jw
    103:   PetscErrorCode ierr
    104:   KSP ksp
    105:   DM dm

    107:   PetscCallA(PetscInitialize(ierr))
    108:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
    109:   PetscCallA(DMDACreate2D(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 3_PETSC_INT_KIND, 3_PETSC_INT_KIND, PETSC_DECIDE, PETSC_DECIDE, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
    110:   PetscCallA(DMSetFromOptions(dm, ierr))
    111:   PetscCallA(DMSetUp(dm, ierr))
    112:   PetscCallA(KSPSetDM(ksp, dm, ierr))
    113:   PetscCallA(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, 0, ierr))
    114:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, 0, ierr))
    115:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, 0, ierr))
    116:   PetscCallA(DMDAGetCorners(dm, is, js, PETSC_NULL_INTEGER, iw, jw, PETSC_NULL_INTEGER, ierr))
    117:   PetscCallA(KSPSetFromOptions(ksp, ierr))
    118:   PetscCallA(KSPSetUp(ksp, ierr))
    119:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
    120:   PetscCallA(KSPDestroy(ksp, ierr))
    121:   PetscCallA(DMDestroy(dm, ierr))
    122:   PetscCallA(PetscFinalize(ierr))
    123: end

    125: !/*TEST
    126: !
    127: !   test:
    128: !      nsize: 4
    129: !      args: -ksp_monitor_short -da_refine 5 -pc_type mg -pc_mg_levels 5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -mg_levels_pc_type jacobi -ksp_pc_side right
    130: !
    131: !TEST*/