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 ierr
14: KSP ksp
15: PetscInt ctx(*)
16: Vec b
17: PetscScalar h
19: h = 0.0
20: PetscCall(VecSet(b, h, ierr))
21: end subroutine
23: subroutine ComputeRHS(ksp, b, dummy, ierr)
25: PetscErrorCode ierr
26: KSP ksp
27: Vec b
28: integer dummy(*)
29: PetscScalar h, Hx, Hy
30: PetscInt mx, my
31: DM dm
33: PetscCall(KSPGetDM(ksp, dm, ierr))
34: 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))
36: Hx = 1.0/real(mx - 1)
37: Hy = 1.0/real(my - 1)
38: h = Hx*Hy
39: PetscCall(VecSet(b, h, ierr))
40: end subroutine
42: subroutine ComputeMatrix(ksp, A, B, dummy, ierr)
43: PetscErrorCode ierr
44: KSP ksp
45: Mat A, B
46: integer dummy(*)
47: DM dm
49: PetscInt i, j, mx, my, xm
50: PetscInt ym, xs, ys, i1, i5
51: PetscScalar v(5), Hx, Hy
52: PetscScalar HxdHy, HydHx
53: MatStencil row(1), col(5)
55: i1 = 1
56: i5 = 5
57: PetscCall(KSPGetDM(ksp, dm, ierr))
58: 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))
60: Hx = 1.0/real(mx - 1)
61: Hy = 1.0/real(my - 1)
62: HxdHy = Hx/Hy
63: HydHx = Hy/Hx
64: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
65: do j = ys, ys + ym - 1
66: do i = xs, xs + xm - 1
67: row(1)%i = i
68: row(1)%j = j
69: if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
70: v(1) = 2.0*(HxdHy + HydHx)
71: PetscCall(MatSetValuesStencil(B, i1, row, i1, row, v, INSERT_VALUES, ierr))
72: else
73: v(1) = -HxdHy
74: col(1)%i = i
75: col(1)%j = j - 1
76: v(2) = -HydHx
77: col(2)%i = i - 1
78: col(2)%j = j
79: v(3) = 2.0*(HxdHy + HydHx)
80: col(3)%i = i
81: col(3)%j = j
82: v(4) = -HydHx
83: col(4)%i = i + 1
84: col(4)%j = j
85: v(5) = -HxdHy
86: col(5)%i = i
87: col(5)%j = j + 1
88: PetscCall(MatSetValuesStencil(B, i1, row, i5, col, v, INSERT_VALUES, ierr))
89: end if
90: end do
91: end do
92: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
93: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
94: if (A /= B) then
95: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
96: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
97: end if
98: end subroutine
99: end module ex45fmodule
101: program main
102: use petscdmda
103: use petscksp
104: use ex45fmodule
105: implicit none
107: PetscInt is, js, iw, jw
108: PetscInt one, three
109: PetscErrorCode ierr
110: KSP ksp
111: DM dm
113: one = 1
114: three = 3
116: PetscCallA(PetscInitialize(ierr))
117: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
118: PetscCallA(DMDACreate2D(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, three, three, PETSC_DECIDE, PETSC_DECIDE, one, one, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
119: PetscCallA(DMSetFromOptions(dm, ierr))
120: PetscCallA(DMSetUp(dm, ierr))
121: PetscCallA(KSPSetDM(ksp, dm, ierr))
122: PetscCallA(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, 0, ierr))
123: PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, 0, ierr))
124: PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, 0, ierr))
125: PetscCallA(DMDAGetCorners(dm, is, js, PETSC_NULL_INTEGER, iw, jw, PETSC_NULL_INTEGER, ierr))
126: PetscCallA(KSPSetFromOptions(ksp, ierr))
127: PetscCallA(KSPSetUp(ksp, ierr))
128: PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
129: PetscCallA(KSPDestroy(ksp, ierr))
130: PetscCallA(DMDestroy(dm, ierr))
131: PetscCallA(PetscFinalize(ierr))
132: end
134: !/*TEST
135: !
136: ! test:
137: ! nsize: 4
138: ! 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
139: !
140: !TEST*/