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 12 program main 13 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,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,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 46 end 47 48 subroutine 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)) 65 return 66 end 67 68 subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr) 69 use petscksp 70 implicit none 71 72 Mat jac,JJ 73 PetscErrorCode ierr 74 KSP ksp 75 DM da 76 PetscInt i,j,k,mx,my,mz,xm 77 PetscInt ym,zm,xs,ys,zs,i1,i7 78 PetscScalar v(7),Hx,Hy,Hz 79 PetscScalar HxHydHz,HyHzdHx 80 PetscScalar HxHzdHy 81 MatStencil row(4),col(4,7) 82 PetscInt ctx 83 i1 = 1 84 i7 = 7 85 PetscCall(KSPGetDM(ksp,da,ierr)) 86 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)) 87 88 Hx = 1.0 / real(mx-1) 89 Hy = 1.0 / real(my-1) 90 Hz = 1.0 / real(mz-1) 91 HxHydHz = Hx*Hy/Hz 92 HxHzdHy = Hx*Hz/Hy 93 HyHzdHx = Hy*Hz/Hx 94 PetscCall(DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)) 95 96 do 10,k=zs,zs+zm-1 97 do 20,j=ys,ys+ym-1 98 do 30,i=xs,xs+xm-1 99 row(MatStencil_i) = i 100 row(MatStencil_j) = j 101 row(MatStencil_k) = k 102 if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 .or. k.eq.mz-1) then 103 v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx) 104 PetscCall(MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,ierr)) 105 else 106 v(1) = -HxHydHz 107 col(MatStencil_i,1) = i 108 col(MatStencil_j,1) = j 109 col(MatStencil_k,1) = k-1 110 v(2) = -HxHzdHy 111 col(MatStencil_i,2) = i 112 col(MatStencil_j,2) = j-1 113 col(MatStencil_k,2) = k 114 v(3) = -HyHzdHx 115 col(MatStencil_i,3) = i-1 116 col(MatStencil_j,3) = j 117 col(MatStencil_k,3) = k 118 v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx) 119 col(MatStencil_i,4) = i 120 col(MatStencil_j,4) = j 121 col(MatStencil_k,4) = k 122 v(5) = -HyHzdHx 123 col(MatStencil_i,5) = i+1 124 col(MatStencil_j,5) = j 125 col(MatStencil_k,5) = k 126 v(6) = -HxHzdHy 127 col(MatStencil_i,6) = i 128 col(MatStencil_j,6) = j+1 129 col(MatStencil_k,6) = k 130 v(7) = -HxHydHz 131 col(MatStencil_i,7) = i 132 col(MatStencil_j,7) = j 133 col(MatStencil_k,7) = k+1 134 PetscCall(MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES,ierr)) 135 endif 136 30 continue 137 20 continue 138 10 continue 139 140 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 141 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 142 return 143 end 144 145!/*TEST 146! 147! test: 148! 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 149! requires: !single 150! output_file: output/ex22_1.out 151! 152!TEST*/ 153