Actual source code: ex54f.F90

      1: !
      2: !   Description: Solve Ax=b.  A comes from an anisotropic 2D thermal problem with Q1 FEM on domain (-1,1)^2.
      3: !       Material conductivity given by tensor:
      4: !
      5: !       D = | 1 0       |
      6: !           | 0 epsilon |
      7: !
      8: !    rotated by angle 'theta' (-theta <90> in degrees) with anisotropic parameter 'epsilon' (-epsilon <0.0>).
      9: !    Blob right-hand side centered at C (-blob_center C(1),C(2) <0,0>)
     10: !    Dirichlet BCs on y=-1 face.
     11: !
     12: !    -out_matlab will generate binary files for A,x,b and a ex54f.m file that reads them and plots them in matlab.
     13: !
     14: !    User can change anisotropic shape with function ex54_psi().  Negative theta will switch to a circular anisotropy.
     15: !

     17: ! -----------------------------------------------------------------------
     18: #include <petsc/finclude/petscksp.h>
     19: module ex54fmodule
     20:   use petscksp

     22:   implicit none
     23:   PetscReal theta

     25: contains
     26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     27: !     thfx2d - compute material tensor
     28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     29: !     Compute thermal gradient and flux

     31:   pure subroutine thfx2d(ev, xl, shp, dd, ndm, nel, dir)
     32:     PetscInt, intent(in) :: nel, ndm
     33:     PetscReal, intent(out) :: dd(2, 2)
     34:     PetscReal, intent(in) :: ev(2), xl(ndm, nel), shp(3, *)
     35:     procedure(ex54_psi) :: dir
     36:     PetscReal xx, yy, psi, cs, sn, c2, s2

     38:     xx = sum(shp(3, 1:nel)*xl(1, 1:nel))
     39:     yy = sum(shp(3, 1:nel)*xl(2, 1:nel))
     40:     psi = dir(xx, yy)
     41: !   Compute thermal flux
     42:     cs = cos(psi)
     43:     sn = sin(psi)
     44:     c2 = cs*cs
     45:     s2 = sn*sn
     46:     cs = cs*sn

     48:     dd(1, 1) = c2*ev(1) + s2*ev(2)
     49:     dd(2, 2) = s2*ev(1) + c2*ev(2)
     50:     dd(1, 2) = cs*(ev(1) - ev(2))
     51:     dd(2, 1) = dd(1, 2)

     53: !   flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2)
     54: !   flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2)
     55:   end subroutine thfx2d

     57: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     58: !     shp2dquad - shape functions - compute derivatives w/r natural coords.
     59: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60:   pure subroutine shp2dquad(s, t, xl, shp, xsj, ndm)
     61: !-----[--.----+----.----+----.-----------------------------------------]
     62: !      Purpose: Shape function routine for 4-node isoparametric quads
     63: !
     64: !      Inputs:
     65: !         s,t       - Natural coordinates of point
     66: !         xl(ndm,*) - Nodal coordinates for element
     67: !         ndm       - Spatial dimension of mesh

     69: !      Outputs:
     70: !         shp(3,*)  - Shape functions and derivatives at point
     71: !                     shp(1,i) = dN_i/dx  or dN_i/dxi_1
     72: !                     shp(2,i) = dN_i/dy  or dN_i/dxi_2
     73: !                     shp(3,i) = N_i
     74: !         xsj       - Jacobian determinant at point
     75: !-----[--.----+----.----+----.-----------------------------------------]
     76:     PetscInt, intent(in) :: ndm
     77:     PetscReal, intent(out) :: shp(3, 4), xsj
     78:     PetscReal, intent(in) :: xl(ndm, 4), s, t
     79:     PetscReal xo, xs, xt, yo, ys, yt, xsm, xsp, xtm
     80:     PetscReal xtp, ysm, ysp, ytm, ytp
     81:     PetscReal xsj1, sh, th, sp, tp, sm, tm

     83: !   Set up interpolations

     85:     sh = 0.5*s
     86:     th = 0.5*t
     87:     sp = 0.5 + sh
     88:     tp = 0.5 + th
     89:     sm = 0.5 - sh
     90:     tm = 0.5 - th
     91:     shp(3, 1) = sm*tm
     92:     shp(3, 2) = sp*tm
     93:     shp(3, 3) = sp*tp
     94:     shp(3, 4) = sm*tp

     96: !   Set up natural coordinate functions (times 4)

     98:     xo = xl(1, 1) - xl(1, 2) + xl(1, 3) - xl(1, 4)
     99:     xs = -xl(1, 1) + xl(1, 2) + xl(1, 3) - xl(1, 4) + xo*t
    100:     xt = -xl(1, 1) - xl(1, 2) + xl(1, 3) + xl(1, 4) + xo*s
    101:     yo = xl(2, 1) - xl(2, 2) + xl(2, 3) - xl(2, 4)
    102:     ys = -xl(2, 1) + xl(2, 2) + xl(2, 3) - xl(2, 4) + yo*t
    103:     yt = -xl(2, 1) - xl(2, 2) + xl(2, 3) + xl(2, 4) + yo*s

    105: !   Compute jacobian (times 16)

    107:     xsj1 = xs*yt - xt*ys

    109: !   Divide jacobian by 16 (multiply by .0625)

    111:     xsj = 0.0625*xsj1
    112:     if (xsj1 == 0.0) then
    113:       xsj1 = 1.0
    114:     else
    115:       xsj1 = 1.0/xsj1
    116:     end if

    118: !   Divide functions by jacobian

    120:     xs = (xs + xs)*xsj1
    121:     xt = (xt + xt)*xsj1
    122:     ys = (ys + ys)*xsj1
    123:     yt = (yt + yt)*xsj1

    125: !   Multiply by interpolations

    127:     ytm = yt*tm
    128:     ysm = ys*sm
    129:     ytp = yt*tp
    130:     ysp = ys*sp
    131:     xtm = xt*tm
    132:     xsm = xs*sm
    133:     xtp = xt*tp
    134:     xsp = xs*sp

    136: !   Compute shape functions

    138:     shp(1, 1) = -ytm + ysm
    139:     shp(1, 2) = ytm + ysp
    140:     shp(1, 3) = ytp - ysp
    141:     shp(1, 4) = -ytp - ysm
    142:     shp(2, 1) = xtm - xsm
    143:     shp(2, 2) = -xtm - xsp
    144:     shp(2, 3) = -xtp + xsp
    145:     shp(2, 4) = xtp + xsm

    147:   end subroutine shp2dquad

    149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    150: ! int2d
    151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    152:   subroutine int2d(l, sg)
    153: !-----[--.----+----.----+----.-----------------------------------------]
    154: !     Purpose: Form Gauss points and weights for two dimensions

    156: !     Inputs:
    157: !     l       - Number of points/direction

    159: !     Outputs:
    160: !     sg(3,*) - Array of points and weights
    161: !-----[--.----+----.----+----.-----------------------------------------]
    162:     PetscInt l, i
    163:     PetscReal g, sg(3, *)
    164:     PetscInt, parameter, dimension(9) :: &
    165:       lr = [-1, 1, 1, -1, 0, 1, 0, -1, 0], &
    166:       lz = [-1, -1, 1, 1, -1, 0, 1, 0, 0]
    167:     PetscReal, parameter :: third = 1.0_PETSC_REAL_KIND/3.0_PETSC_REAL_KIND

    169: !   2x2 integration
    170:     g = sqrt(third)
    171:     do i = 1, 4
    172:       sg(1, i) = g*lr(i)
    173:       sg(2, i) = g*lz(i)
    174:       sg(3, i) = 1.0
    175:     end do

    177:   end

    179: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    180: ! ex54_psi - anisotropic material direction
    181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    182:   PetscReal pure function ex54_psi(x, y)
    183:     PetscReal, intent(in) :: x, y
    184:     ex54_psi = theta
    185:     if (theta < 0.) then     ! circular
    186:       if (y == 0.0) then
    187:         ex54_psi = 2.0*atan(1.0)
    188:       else
    189:         ex54_psi = atan(-x/y)
    190:       end if
    191:     end if
    192:   end
    193: end module ex54fmodule

    195: program main
    196:   use petscksp
    197:   use ex54fmodule
    198:   implicit none

    200:   Vec xvec, bvec, uvec
    201:   Mat Amat
    202:   KSP ksp
    203:   PetscErrorCode ierr
    204:   PetscViewer viewer
    205:   PetscInt qj, qi, ne, M, Istart, Iend, geq
    206:   PetscInt ki, kj, ll, j1, i1
    207:   PetscInt, parameter :: ndf = 1
    208:   PetscInt, parameter :: nel = 4 ! nodes per element (quad)
    209:   PetscInt :: idx(4)
    210:   PetscBool flg, out_matlab
    211:   PetscMPIInt size, rank
    212:   PetscScalar :: ss(4, 4), val
    213:   PetscReal :: shp(3, 9), sg(3, 9)
    214:   PetscReal :: thk = 1.0              ! thickness
    215:   PetscReal :: eps, h, x, y, xsj, a1, a2
    216:   PetscReal :: coord(2, 4), dd(2, 2), ev(3), blb(2)
    217:   PetscReal, parameter :: rad2deg = 57.2957795_PETSC_REAL_KIND

    219: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    220: !                 Beginning of program
    221: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    222:   PetscCallA(PetscInitialize(ierr))
    223:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
    224:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    225: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    226: !                 set parameters
    227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    228:   ne = 9
    229:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-ne', ne, flg, ierr))
    230:   h = 2.0/real(ne)
    231:   M = (ne + 1)**2
    232:   theta = 90.0
    233: ! theta is input in degrees
    234:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-theta', theta, flg, ierr))
    235:   theta = theta/rad2deg
    236:   eps = 1.0
    237:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-epsilon', eps, flg, ierr))
    238:   ki = 2
    239:   PetscCallA(PetscOptionsGetRealArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-blob_center', blb, ki, flg, ierr))
    240:   if (.not. flg) then
    241:     blb = 0.0
    242:   else if (ki /= 2) then
    243:     print *, 'error: ', ki, ' arguments read for -blob_center.  Needs to be two.'
    244:   end if
    245:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-out_matlab', out_matlab, flg, ierr))
    246:   if (.not. flg) out_matlab = PETSC_FALSE

    248:   ev(1) = 1.0
    249:   ev(2) = eps*ev(1)
    250: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    251: !    Compute the matrix and right-hand-side vector that define
    252: !    the linear system, Ax = b.
    253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    254: ! Create matrix.  When using MatCreate(), the matrix format can
    255: ! be specified at runtime.
    256:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr))
    257:   PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, M, M, ierr))
    258:   PetscCallA(MatSetType(Amat, MATAIJ, ierr))
    259:   PetscCallA(MatSetOption(Amat, MAT_SPD, PETSC_TRUE, ierr))
    260:   PetscCallA(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE, ierr))
    261:   if (size == 1) then
    262:     PetscCallA(MatSetType(Amat, MATAIJ, ierr))
    263:   else
    264:     PetscCallA(MatSetType(Amat, MATMPIAIJ, ierr))
    265:   end if
    266:   PetscCallA(MatMPIAIJSetPreallocation(Amat, 9_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 6_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    267:   PetscCallA(MatSetFromOptions(Amat, ierr))
    268:   PetscCallA(MatSetUp(Amat, ierr))
    269:   PetscCallA(MatGetOwnershipRange(Amat, Istart, Iend, ierr))
    270: !  Create vectors.  Note that we form 1 vector from scratch and
    271: !  then duplicate as needed.
    272:   PetscCallA(MatCreateVecs(Amat, PETSC_NULL_VEC, xvec, ierr))
    273:   PetscCallA(VecSetFromOptions(xvec, ierr))
    274:   PetscCallA(VecDuplicate(xvec, bvec, ierr))
    275:   PetscCallA(VecDuplicate(xvec, uvec, ierr))
    276: ! Assemble matrix.
    277: !  - Note that MatSetValues() uses 0-based row and column numbers
    278: !    in Fortran as well as in C (as set here in the array "col").
    279:   call int2d(2_PETSC_INT_KIND, sg)
    280:   do geq = Istart, Iend - 1, 1
    281:     qj = geq/(ne + 1); qi = mod(geq, (ne + 1))
    282:     x = h*qi - 1.0; y = h*qj - 1.0 ! lower left corner (-1,-1)
    283:     if (qi < ne .and. qj < ne) then
    284:       coord(1, 1) = x; coord(2, 1) = y
    285:       coord(1, 2) = x + h; coord(2, 2) = y
    286:       coord(1, 3) = x + h; coord(2, 3) = y + h
    287:       coord(1, 4) = x; coord(2, 4) = y + h
    288: ! form stiff
    289:       ss = 0.0
    290:       do ll = 1, 4
    291:         call shp2dquad(sg(1, ll), sg(2, ll), coord, shp, xsj, 2_PETSC_INT_KIND)
    292:         xsj = xsj*sg(3, ll)*thk
    293:         call thfx2d(ev, coord, shp, dd, 2_PETSC_INT_KIND, 4_PETSC_INT_KIND, ex54_psi)
    294:         j1 = 1
    295:         do kj = 1, nel
    296:           a1 = (dd(1, 1)*shp(1, kj) + dd(1, 2)*shp(2, kj))*xsj
    297:           a2 = (dd(2, 1)*shp(1, kj) + dd(2, 2)*shp(2, kj))*xsj
    298: !     Compute residual
    299: !         p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2)
    300: !     Compute tangent
    301:           i1 = 1
    302:           do ki = 1, nel
    303:             ss(i1, j1) = ss(i1, j1) + a1*shp(1, ki) + a2*shp(2, ki)
    304:             i1 = i1 + ndf
    305:           end do
    306:           j1 = j1 + ndf
    307:         end do
    308:       end do

    310:       idx = [geq, geq + 1, geq + (ne + 1) + 1, geq + (ne + 1)]
    311:       if (qj > 0) then
    312:         PetscCallA(MatSetValues(Amat, 4_PETSC_INT_KIND, idx, 4_PETSC_INT_KIND, idx, reshape(ss, [4_PETSC_INT_KIND**2]), ADD_VALUES, ierr))
    313:       else                !     a BC
    314:         do ki = 1, 4, 1
    315:           do kj = 1, 4, 1
    316:             if (ki < 3 .or. kj < 3) then
    317:               if (ki == kj) then
    318:                 ss(ki, kj) = .1*ss(ki, kj)
    319:               else
    320:                 ss(ki, kj) = 0.0
    321:               end if
    322:             end if
    323:           end do
    324:         end do
    325:         PetscCallA(MatSetValues(Amat, 4_PETSC_INT_KIND, idx, 4_PETSC_INT_KIND, idx, reshape(ss, [4_PETSC_INT_KIND**2]), ADD_VALUES, ierr))
    326:       end if               ! BC
    327:     end if                  ! add element
    328:     if (qj > 0) then      ! set rhs
    329:       val = h*h*exp(-100*((x + h/2) - blb(1))**2)*exp(-100*((y + h/2) - blb(2))**2)
    330:       PetscCallA(VecSetValues(bvec, 1_PETSC_INT_KIND, [geq], [val], INSERT_VALUES, ierr))
    331:     end if
    332:   end do
    333:   PetscCallA(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY, ierr))
    334:   PetscCallA(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY, ierr))
    335:   PetscCallA(VecAssemblyBegin(bvec, ierr))
    336:   PetscCallA(VecAssemblyEnd(bvec, ierr))

    338: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    339: !          Create the linear solver and set various options
    340: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    342: !  Create linear solver context

    344:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))

    346: !  Set operators. Here the matrix that defines the linear system
    347: !  also serves as the matrix from which the preconditioner is constructed.

    349:   PetscCallA(KSPSetOperators(ksp, Amat, Amat, ierr))

    351: !  Set runtime options, e.g.,
    352: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    353: !  These options will override those specified above as long as
    354: !  KSPSetFromOptions() is called _after_ any other customization
    355: !  routines.

    357:   PetscCallA(KSPSetFromOptions(ksp, ierr))

    359: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    360: !                      Solve the linear system
    361: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    363:   PetscCallA(KSPSolve(ksp, bvec, xvec, ierr))

    365: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    366: !                      output
    367: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    368:   if (out_matlab) then
    369:     PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Amat', FILE_MODE_WRITE, viewer, ierr))
    370:     PetscCallA(MatView(Amat, viewer, ierr))
    371:     PetscCallA(PetscViewerDestroy(viewer, ierr))

    373:     PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Bvec', FILE_MODE_WRITE, viewer, ierr))
    374:     PetscCallA(VecView(bvec, viewer, ierr))
    375:     PetscCallA(PetscViewerDestroy(viewer, ierr))

    377:     PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Xvec', FILE_MODE_WRITE, viewer, ierr))
    378:     PetscCallA(VecView(xvec, viewer, ierr))
    379:     PetscCallA(PetscViewerDestroy(viewer, ierr))

    381:     PetscCallA(MatMult(Amat, xvec, uvec, ierr))
    382:     val = -1.0
    383:     PetscCallA(VecAXPY(uvec, val, bvec, ierr))
    384:     PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Rvec', FILE_MODE_WRITE, viewer, ierr))
    385:     PetscCallA(VecView(uvec, viewer, ierr))
    386:     PetscCallA(PetscViewerDestroy(viewer, ierr))

    388:     if (rank == 0) then
    389:       open (1, file='ex54f.m', FORM='formatted')
    390:       write (1, *) 'A = PetscBinaryRead(''Amat'');'
    391:       write (1, *) '[m n] = size(A);'
    392:       write (1, *) 'mm = sqrt(m);'
    393:       write (1, *) 'b = PetscBinaryRead(''Bvec'');'
    394:       write (1, *) 'x = PetscBinaryRead(''Xvec'');'
    395:       write (1, *) 'r = PetscBinaryRead(''Rvec'');'
    396:       write (1, *) 'bb = reshape(b,mm,mm);'
    397:       write (1, *) 'xx = reshape(x,mm,mm);'
    398:       write (1, *) 'rr = reshape(r,mm,mm)'
    399: !     write (1,*) 'imagesc(bb')'
    400: !     write (1,*) 'title('RHS'),'
    401:       write (1, *) 'figure,'
    402:       write (1, *) 'imagesc(xx'')'
    403:       write (1, 2002) eps, theta*57.2957795
    404:       write (1, *) 'figure,'
    405:       write (1, *) 'imagesc(rr'')'
    406:       write (1, *) 'title(''Residual''),'
    407:       close (1)
    408:     end if
    409:   end if
    410: 2002 format('title(''Solution: esp='',d9.3,'', theta='',g8.3,''),')
    411: !  Free work space.  All PETSc objects should be destroyed when they
    412: !  are no longer needed.

    414:   PetscCallA(VecDestroy(xvec, ierr))
    415:   PetscCallA(VecDestroy(bvec, ierr))
    416:   PetscCallA(VecDestroy(uvec, ierr))
    417:   PetscCallA(MatDestroy(Amat, ierr))
    418:   PetscCallA(KSPDestroy(ksp, ierr))
    419:   PetscCallA(PetscFinalize(ierr))
    420: end

    422: !
    423: !/*TEST
    424: !
    425: !  testset:
    426: !   nsize: 4
    427: !   args: -ne 39 -theta 30.0 -epsilon 1.e-1 -blob_center 0.,0. -ksp_type cg -pc_type gamg -pc_gamg_type agg -ksp_rtol 1e-4 -pc_gamg_aggressive_square_graph false -ksp_norm_type unpreconditioned
    428: !   requires: !single
    429: !   test:
    430: !      suffix: misk
    431: !      args: -pc_gamg_mat_coarsen_type misk -pc_gamg_aggressive_coarsening 0 -ksp_monitor_short
    432: !   test:
    433: !      suffix: mis
    434: !      args: -pc_gamg_mat_coarsen_type mis -ksp_monitor_short
    435: !   test:
    436: !      suffix: hem
    437: !      args: -pc_gamg_mat_coarsen_type hem -ksp_converged_reason
    438: !      filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 1[2-3]/Linear solve converged due to CONVERGED_RTOL iterations 11/g"
    439: !
    440: !TEST*/