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