1: !
  2: !   Laplacian in 3D. Modeled by the partial differential equation
  3: !
  4: !   Laplacian u = 1,0 < x,y,z < 1,
  5: !
  6: !   with boundary conditions
  7: !
  8: !   u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
  9: !
 10: !   This uses multigrid to solve the linear system
 11: #include <petsc/finclude/petscdmda.h>
 12: #include <petsc/finclude/petscksp.h>
 13: module ex22fmodule
 14:   use petscksp
 15:   use petscdmda
 16:   implicit none

 18: contains
 19:   subroutine ComputeRHS(ksp, b, ctx, ierr)
 20:     PetscErrorCode, intent(out) :: ierr
 21:     PetscInt mx, my, mz
 22:     PetscScalar h
 23:     Vec b
 24:     KSP ksp
 25:     DM da
 26:     PetscInt ctx

 28:     PetscCall(KSPGetDM(ksp, da, ierr))
 29:     PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, my, mz, 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))
 30:     h = 1.0/real((mx - 1)*(my - 1)*(mz - 1))

 32:     PetscCall(VecSet(b, h, ierr))
 33:   end

 35:   subroutine ComputeMatrix(ksp, JJ, jac, ctx, ierr)
 36:     Mat jac, JJ
 37:     PetscErrorCode ierr
 38:     KSP ksp
 39:     DM da
 40:     PetscInt i, j, k, mx, my, mz, xm
 41:     PetscInt ym, zm, xs, ys, zs
 42:     PetscScalar v(7), Hx, Hy, Hz
 43:     PetscScalar HxHydHz, HyHzdHx
 44:     PetscScalar HxHzdHy
 45:     MatStencil row(1), col(7)
 46:     PetscInt ctx
 47:     PetscCall(KSPGetDM(ksp, da, ierr))
 48:     PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, my, mz, 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))

 50:     Hx = 1.0/real(mx - 1)
 51:     Hy = 1.0/real(my - 1)
 52:     Hz = 1.0/real(mz - 1)
 53:     HxHydHz = Hx*Hy/Hz
 54:     HxHzdHy = Hx*Hz/Hy
 55:     HyHzdHx = Hy*Hz/Hx
 56:     PetscCall(DMDAGetCorners(da, xs, ys, zs, xm, ym, zm, ierr))

 58:     do k = zs, zs + zm - 1
 59:       do j = ys, ys + ym - 1
 60:         do i = xs, xs + xm - 1
 61:           row(1)%i = i
 62:           row(1)%j = j
 63:           row(1)%k = k
 64:           if (i == 0 .or. j == 0 .or. k == 0 .or. i == mx - 1 .or. j == my - 1 .or. k == mz - 1) then
 65:             v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
 66:             PetscCall(MatSetValuesStencil(jac, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, row, v, INSERT_VALUES, ierr))
 67:           else
 68:             v(1) = -HxHydHz
 69:             col(1)%i = i
 70:             col(1)%j = j
 71:             col(1)%k = k - 1
 72:             v(2) = -HxHzdHy
 73:             col(2)%i = i
 74:             col(2)%j = j - 1
 75:             col(2)%k = k
 76:             v(3) = -HyHzdHx
 77:             col(3)%i = i - 1
 78:             col(3)%j = j
 79:             col(3)%k = k
 80:             v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
 81:             col(4)%i = i
 82:             col(4)%j = j
 83:             col(4)%k = k
 84:             v(5) = -HyHzdHx
 85:             col(5)%i = i + 1
 86:             col(5)%j = j
 87:             col(5)%k = k
 88:             v(6) = -HxHzdHy
 89:             col(6)%i = i
 90:             col(6)%j = j + 1
 91:             col(6)%k = k
 92:             v(7) = -HxHydHz
 93:             col(7)%i = i
 94:             col(7)%j = j
 95:             col(7)%k = k + 1
 96:             PetscCall(MatSetValuesStencil(jac, 1_PETSC_INT_KIND, row, 7_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
 97:           end if
 98:         end do
 99:       end do
100:     end do

102:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
103:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
104:   end
105: end module ex22fmodule

107: program main
108:   use ex22fmodule
109:   implicit none

111:   PetscErrorCode ierr
112:   DM da
113:   KSP ksp
114:   Vec x
115:   PetscInt ctx

117:   PetscCallA(PetscInitialize(ierr))

119:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
120:   PetscCallA(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 3_PETSC_INT_KIND, 3_PETSC_INT_KIND, 3_PETSC_INT_KIND, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr))
121:   PetscCallA(DMSetFromOptions(da, ierr))
122:   PetscCallA(DMSetUp(da, ierr))
123:   PetscCallA(KSPSetDM(ksp, da, ierr))
124:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, ctx, ierr))
125:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, ctx, ierr))

127:   PetscCallA(KSPSetFromOptions(ksp, ierr))
128:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
129:   PetscCallA(KSPGetSolution(ksp, x, ierr))
130:   PetscCallA(KSPDestroy(ksp, ierr))
131:   PetscCallA(DMDestroy(da, ierr))
132:   PetscCallA(PetscFinalize(ierr))

134: end

136: !/*TEST
137: !
138: !   test:
139: !      args: -pc_mg_type full -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type preconditioned -pc_type mg -da_refine 2 -ksp_type fgmres
140: !      requires: !single
141: !      output_file: output/ex22_1.out
142: !
143: !TEST*/