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