Actual source code: ex22f.F90

  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

 14: #include <petsc/finclude/petscksp.h>
 15:       use petscdmda
 16:       use petscksp
 17:       implicit none

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

 27:       PetscCallA(PetscInitialize(ierr))

 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))

 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))

 46:       end

 48:       subroutine ComputeRHS(ksp,b,ctx,ierr)
 49:       use petscksp
 50:       implicit none

 52:       PetscErrorCode  ierr
 53:       PetscInt mx,my,mz
 54:       PetscScalar  h
 55:       Vec          b
 56:       KSP          ksp
 57:       DM           da
 58:       PetscInt     ctx

 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))

 64:       PetscCall(VecSet(b,h,ierr))
 65:       return
 66:       end

 68:       subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr)
 69:       use petscksp
 70:       implicit none

 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))

 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))

 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

140:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
141:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
142:       return
143:       end

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*/