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