Actual source code: ex22f.F

petsc-3.8.3 2017-12-09
Report Typos and Errors
  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

 12:       program main
 13:  #include <petsc/finclude/petscksp.h>
 14:       use petscdmda
 15:       use petscksp
 16:       implicit none

 18:       PetscErrorCode   ierr
 19:       DM               da
 20:       KSP              ksp
 21:       Vec              x
 22:       external         ComputeRHS,ComputeMatrix
 23:       PetscInt i1,i3
 24:       PetscInt         ctx

 26:       call  PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 27:       if (ierr .ne. 0) then
 28:         print*,'Unable to initialize PETSc'
 29:         stop
 30:       endif

 32:       i3 = 3
 33:       i1 = 1
 34:       call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
 35:       call DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,              &
 36:      &    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,                        &
 39:      &                PETSC_NULL_INTEGER,da,ierr)
 40:       call DMSetFromOptions(da,ierr)
 41:       call DMSetUp(da,ierr)
 42:       call KSPSetDM(ksp,da,ierr)
 43:       call KSPSetComputeRHS(ksp,ComputeRHS,ctx,ierr)
 44:       call KSPSetComputeOperators(ksp,ComputeMatrix,ctx,ierr)

 46:       call KSPSetFromOptions(ksp,ierr)
 47:       call KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
 48:       call KSPGetSolution(ksp,x,ierr)
 49:       call KSPDestroy(ksp,ierr)
 50:       call DMDestroy(da,ierr)
 51:       call PetscFinalize(ierr)

 53:       end

 55:       subroutine ComputeRHS(ksp,b,ctx,ierr)
 56:       use petscksp
 57:       implicit none

 59:       PetscErrorCode  ierr
 60:       PetscInt mx,my,mz
 61:       PetscScalar  h
 62:       Vec          b
 63:       KSP          ksp
 64:       DM           da
 65:       PetscInt     ctx

 67:       call KSPGetDM(ksp,da,ierr)
 68:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                        &
 69:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 70:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 71:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 72:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 73:      &               PETSC_NULL_INTEGER,ierr)
 74:       h    = 1.0/real((mx-1)*(my-1)*(mz-1))

 76:       call VecSet(b,h,ierr)
 77:       return
 78:       end


 81:       subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr)
 82:       use petscksp
 83:       implicit none

 85:       Mat          jac,JJ
 86:       PetscErrorCode    ierr
 87:       KSP          ksp
 88:       DM           da
 89:       PetscInt    i,j,k,mx,my,mz,xm
 90:       PetscInt    ym,zm,xs,ys,zs,i1,i7
 91:       PetscScalar  v(7),Hx,Hy,Hz
 92:       PetscScalar  HxHydHz,HyHzdHx
 93:       PetscScalar  HxHzdHy
 94:       MatStencil   row(4),col(4,7)
 95:       PetscInt     ctx
 96:       i1 = 1
 97:       i7 = 7
 98:       call KSPGetDM(ksp,da,ierr)
 99:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                       &
100:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
101:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
102:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
103:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
104:      &               PETSC_NULL_INTEGER,ierr)

106:       Hx = 1.0 / real(mx-1)
107:       Hy = 1.0 / real(my-1)
108:       Hz = 1.0 / real(mz-1)
109:       HxHydHz = Hx*Hy/Hz
110:       HxHzdHy = Hx*Hz/Hy
111:       HyHzdHx = Hy*Hz/Hx
112:       call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)

114:       do 10,k=zs,zs+zm-1
115:         do 20,j=ys,ys+ym-1
116:           do 30,i=xs,xs+xm-1
117:           row(MatStencil_i) = i
118:           row(MatStencil_j) = j
119:           row(MatStencil_k) = k
120:           if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or.         &
121:      &         j.eq.my-1 .or. k.eq.mz-1) then
122:             v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
123:             call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,    &
124:      &                               ierr)
125:           else
126:             v(1) = -HxHydHz
127:              col(MatStencil_i,1) = i
128:              col(MatStencil_j,1) = j
129:              col(MatStencil_k,1) = k-1
130:             v(2) = -HxHzdHy
131:              col(MatStencil_i,2) = i
132:              col(MatStencil_j,2) = j-1
133:              col(MatStencil_k,2) = k
134:             v(3) = -HyHzdHx
135:              col(MatStencil_i,3) = i-1
136:              col(MatStencil_j,3) = j
137:              col(MatStencil_k,3) = k
138:             v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
139:              col(MatStencil_i,4) = i
140:              col(MatStencil_j,4) = j
141:              col(MatStencil_k,4) = k
142:             v(5) = -HyHzdHx
143:              col(MatStencil_i,5) = i+1
144:              col(MatStencil_j,5) = j
145:              col(MatStencil_k,5) = k
146:             v(6) = -HxHzdHy
147:              col(MatStencil_i,6) = i
148:              col(MatStencil_j,6) = j+1
149:              col(MatStencil_k,6) = k
150:             v(7) = -HxHydHz
151:              col(MatStencil_i,7) = i
152:              col(MatStencil_j,7) = j
153:              col(MatStencil_k,7) = k+1
154:       call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES,               &
155:      &                               ierr)
156:           endif
157:  30       continue
158:  20     continue
159:  10   continue

161:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
162:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
163:       return
164:       end