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