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
     11: #include <petsc/finclude/petscdmda.h>
     12: #include <petsc/finclude/petscksp.h>
     13: module ex22fmodule
     14:   use petscksp
     15:   use petscdmda
     16:   implicit none

     18: contains
     19:   subroutine ComputeRHS(ksp, b, ctx, ierr)
     20:     PetscErrorCode, intent(out) :: ierr
     21:     PetscInt mx, my, mz
     22:     PetscScalar h
     23:     Vec b
     24:     KSP ksp
     25:     DM da
     26:     PetscInt ctx

     28:     PetscCall(KSPGetDM(ksp, da, ierr))
     29:     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_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))
     30:     h = 1.0/real((mx - 1)*(my - 1)*(mz - 1))

     32:     PetscCall(VecSet(b, h, ierr))
     33:   end

     35:   subroutine ComputeMatrix(ksp, JJ, jac, ctx, ierr)
     36:     Mat jac, JJ
     37:     PetscErrorCode ierr
     38:     KSP ksp
     39:     DM da
     40:     PetscInt i, j, k, mx, my, mz, xm
     41:     PetscInt ym, zm, xs, ys, zs
     42:     PetscScalar v(7), Hx, Hy, Hz
     43:     PetscScalar HxHydHz, HyHzdHx
     44:     PetscScalar HxHzdHy
     45:     MatStencil row(1), col(7)
     46:     PetscInt ctx
     47:     PetscCall(KSPGetDM(ksp, da, ierr))
     48:     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_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

     50:     Hx = 1.0/real(mx - 1)
     51:     Hy = 1.0/real(my - 1)
     52:     Hz = 1.0/real(mz - 1)
     53:     HxHydHz = Hx*Hy/Hz
     54:     HxHzdHy = Hx*Hz/Hy
     55:     HyHzdHx = Hy*Hz/Hx
     56:     PetscCall(DMDAGetCorners(da, xs, ys, zs, xm, ym, zm, ierr))

     58:     do k = zs, zs + zm - 1
     59:       do j = ys, ys + ym - 1
     60:         do i = xs, xs + xm - 1
     61:           row(1)%i = i
     62:           row(1)%j = j
     63:           row(1)%k = k
     64:           if (i == 0 .or. j == 0 .or. k == 0 .or. i == mx - 1 .or. j == my - 1 .or. k == mz - 1) then
     65:             v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
     66:             PetscCall(MatSetValuesStencil(jac, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, row, v, INSERT_VALUES, ierr))
     67:           else
     68:             v(1) = -HxHydHz
     69:             col(1)%i = i
     70:             col(1)%j = j
     71:             col(1)%k = k - 1
     72:             v(2) = -HxHzdHy
     73:             col(2)%i = i
     74:             col(2)%j = j - 1
     75:             col(2)%k = k
     76:             v(3) = -HyHzdHx
     77:             col(3)%i = i - 1
     78:             col(3)%j = j
     79:             col(3)%k = k
     80:             v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
     81:             col(4)%i = i
     82:             col(4)%j = j
     83:             col(4)%k = k
     84:             v(5) = -HyHzdHx
     85:             col(5)%i = i + 1
     86:             col(5)%j = j
     87:             col(5)%k = k
     88:             v(6) = -HxHzdHy
     89:             col(6)%i = i
     90:             col(6)%j = j + 1
     91:             col(6)%k = k
     92:             v(7) = -HxHydHz
     93:             col(7)%i = i
     94:             col(7)%j = j
     95:             col(7)%k = k + 1
     96:             PetscCall(MatSetValuesStencil(jac, 1_PETSC_INT_KIND, row, 7_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
     97:           end if
     98:         end do
     99:       end do
    100:     end do

    102:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
    103:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
    104:   end
    105: end module ex22fmodule

    107: program main
    108:   use ex22fmodule
    109:   implicit none

    111:   PetscErrorCode ierr
    112:   DM da
    113:   KSP ksp
    114:   Vec x
    115:   PetscInt ctx

    117:   PetscCallA(PetscInitialize(ierr))

    119:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
    120:   PetscCallA(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 3_PETSC_INT_KIND, 3_PETSC_INT_KIND, 3_PETSC_INT_KIND, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr))
    121:   PetscCallA(DMSetFromOptions(da, ierr))
    122:   PetscCallA(DMSetUp(da, ierr))
    123:   PetscCallA(KSPSetDM(ksp, da, ierr))
    124:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, ctx, ierr))
    125:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, ctx, ierr))

    127:   PetscCallA(KSPSetFromOptions(ksp, ierr))
    128:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
    129:   PetscCallA(KSPGetSolution(ksp, x, ierr))
    130:   PetscCallA(KSPDestroy(ksp, ierr))
    131:   PetscCallA(DMDestroy(da, ierr))
    132:   PetscCallA(PetscFinalize(ierr))

    134: end

    136: !/*TEST
    137: !
    138: !   test:
    139: !      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
    140: !      requires: !single
    141: !      output_file: output/ex22_1.out
    142: !
    143: !TEST*/