Actual source code: ex22f.F

      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:       implicit none

     15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     16: !                    Include files
     17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     18: !
     19: !     petsc.h  - base PETSc routines      petscvec.h - vectors
     20: !     petscsys.h    - system routines          petscmat.h - matrices
     21: !     petscksp.h    - Krylov subspace methods  petscpc.h  - preconditioners

     23:  #include include/finclude/petsc.h
     24:  #include include/finclude/petscvec.h
     25:  #include include/finclude/petscmat.h
     26:  #include include/finclude/petscpc.h
     27:  #include include/finclude/petscksp.h
     28:  #include include/finclude/petscda.h
     29:       DMMG             dmmg
     30:       PetscErrorCode   ierr
     31:       DA               da
     32:       external         ComputeRHS,ComputeJacobian
     33:       Vec              x
     34:       PetscInt i1,i3

     36:       call  PetscInitialize(PETSC_NULL_CHARACTER,ierr)

     38:       i3 = 3
     39:       i1 = 1
     40:       call DMMGCreate(PETSC_COMM_WORLD,i3,PETSC_NULL_INTEGER,dmmg,ierr)
     41:       call DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,            &
     42:      &            i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,i1,i1,               &
     43:      &                PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                        &
     44:      &                PETSC_NULL_INTEGER,da,ierr)
     45:       call DMMGSetDM(dmmg,da,ierr)
     46:       call DADestroy(da,ierr)

     49:       call DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian,ierr)

     51:       call DMMGSolve(dmmg,ierr)

     53:       call DMMGGetx(dmmg,x,ierr)
     54: !      call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)

     56:       call DMMGDestroy(dmmg,ierr)
     57:       call PetscFinalize(ierr)

     59:       end

     61:       subroutine ComputeRHS(dmmg,b,ierr)
     62:       implicit none

     64:  #include include/finclude/petsc.h
     65:  #include include/finclude/petscvec.h
     66:  #include include/finclude/petscda.h

     68:       PetscErrorCode  ierr
     69:       PetscInt mx,my,mz
     70:       PetscScalar  h
     71:       Vec          b
     72:       DMMG         dmmg
     73:       DA           da

     75:       call DMMGGetDA(dmmg,da,ierr)
     76:       call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                        &
     77:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
     78:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
     79:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
     80:      &               PETSC_NULL_INTEGER,ierr)
     81:       h    = 1.d0/((mx-1)*(my-1)*(mz-1))

     83:       call VecSet(b,h,ierr)
     84:       return
     85:       end

     87:
     88:       subroutine ComputeJacobian(dmmg,jac,ierr)
     89:       implicit none

     91:  #include include/finclude/petsc.h
     92:  #include include/finclude/petscvec.h
     93:  #include include/finclude/petscmat.h
     94:  #include include/finclude/petscda.h

     96:       DMMG         dmmg
     97:       Mat          jac
     98:       PetscErrorCode    ierr

    100:       DA           da
    101:       PetscInt    i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,i1,i7
    102:       PetscScalar  v(7),Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy
    103:       MatStencil   row(4),col(4,7)

    105:       i1 = 1
    106:       i7 = 7
    107:       call DMMGGetDA(dmmg,da,ierr)
    108:       call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                       &
    109:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
    110:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
    111:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
    112:      &               PETSC_NULL_INTEGER,ierr)

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

    169:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
    170:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
    171:       return
    172:       end