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> 13module ex22fmodule 14 use petscksp 15 implicit none 16 17contains 18 subroutine ComputeRHS(ksp, b, ctx, ierr) 19 20 PetscErrorCode ierr 21 PetscInt mx, my, mz 22 PetscScalar h 23 Vec b 24 KSP ksp 25 DM da 26 PetscInt ctx 27 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_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)) 30 h = 1.0/real((mx - 1)*(my - 1)*(mz - 1)) 31 32 PetscCall(VecSet(b, h, ierr)) 33 end 34 35 subroutine ComputeMatrix(ksp, JJ, jac, ctx, ierr) 36 37 Mat jac, JJ 38 PetscErrorCode ierr 39 KSP ksp 40 DM da 41 PetscInt i, j, k, mx, my, mz, xm 42 PetscInt ym, zm, xs, ys, zs, i1, i7 43 PetscScalar v(7), Hx, Hy, Hz 44 PetscScalar HxHydHz, HyHzdHx 45 PetscScalar HxHzdHy 46 MatStencil row(1), col(7) 47 PetscInt ctx 48 i1 = 1 49 i7 = 7 50 PetscCall(KSPGetDM(ksp, da, ierr)) 51 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_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)) 52 53 Hx = 1.0/real(mx - 1) 54 Hy = 1.0/real(my - 1) 55 Hz = 1.0/real(mz - 1) 56 HxHydHz = Hx*Hy/Hz 57 HxHzdHy = Hx*Hz/Hy 58 HyHzdHx = Hy*Hz/Hx 59 PetscCall(DMDAGetCorners(da, xs, ys, zs, xm, ym, zm, ierr)) 60 61 do k = zs, zs + zm - 1 62 do j = ys, ys + ym - 1 63 do i = xs, xs + xm - 1 64 row(1)%i = i 65 row(1)%j = j 66 row(1)%k = k 67 if (i == 0 .or. j == 0 .or. k == 0 .or. i == mx - 1 .or. j == my - 1 .or. k == mz - 1) then 68 v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx) 69 PetscCall(MatSetValuesStencil(jac, i1, row, i1, row, v, INSERT_VALUES, ierr)) 70 else 71 v(1) = -HxHydHz 72 col(1)%i = i 73 col(1)%j = j 74 col(1)%k = k - 1 75 v(2) = -HxHzdHy 76 col(2)%i = i 77 col(2)%j = j - 1 78 col(2)%k = k 79 v(3) = -HyHzdHx 80 col(3)%i = i - 1 81 col(3)%j = j 82 col(3)%k = k 83 v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx) 84 col(4)%i = i 85 col(4)%j = j 86 col(4)%k = k 87 v(5) = -HyHzdHx 88 col(5)%i = i + 1 89 col(5)%j = j 90 col(5)%k = k 91 v(6) = -HxHzdHy 92 col(6)%i = i 93 col(6)%j = j + 1 94 col(6)%k = k 95 v(7) = -HxHydHz 96 col(7)%i = i 97 col(7)%j = j 98 col(7)%k = k + 1 99 PetscCall(MatSetValuesStencil(jac, i1, row, i7, col, v, INSERT_VALUES, ierr)) 100 end if 101 end do 102 end do 103 end do 104 105 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 106 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 107 end 108end module ex22fmodule 109 110program main 111 use petscdmda 112 use petscksp 113 use ex22fmodule 114 implicit none 115 116 PetscErrorCode ierr 117 DM da 118 KSP ksp 119 Vec x 120 PetscInt i1, i3 121 PetscInt ctx 122 123 PetscCallA(PetscInitialize(ierr)) 124 125 i3 = 3 126 i1 = 1 127 PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr)) 128 PetscCallA(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, i3, i3, i3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr)) 129 PetscCallA(DMSetFromOptions(da, ierr)) 130 PetscCallA(DMSetUp(da, ierr)) 131 PetscCallA(KSPSetDM(ksp, da, ierr)) 132 PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, ctx, ierr)) 133 PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, ctx, ierr)) 134 135 PetscCallA(KSPSetFromOptions(ksp, ierr)) 136 PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr)) 137 PetscCallA(KSPGetSolution(ksp, x, ierr)) 138 PetscCallA(KSPDestroy(ksp, ierr)) 139 PetscCallA(DMDestroy(da, ierr)) 140 PetscCallA(PetscFinalize(ierr)) 141 142end 143 144!/*TEST 145! 146! test: 147! 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 148! requires: !single 149! output_file: output/ex22_1.out 150! 151!TEST*/ 152