Actual source code: ex73f90t.F90

      1: !
      2: !  Description: Solves a nonlinear system in parallel with SNES.
      3: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
      4: !  domain, using distributed arrays (DMDAs) to partition the parallel grid.
      5: !  The command line options include:
      6: !    -par <parameter>, where <parameter> indicates the nonlinearity of the problem
      7: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
      8: !
      9: !  This system (A) is augmented with constraints:
     10: !
     11: !    A -B   *  phi  =  rho
     12: !   -C  I      lam  = 0
     13: !
     14: !  where I is the identity, A is the "normal" Poisson equation, B is the "distributor" of the
     15: !  total flux (the first block equation is the flux surface averaging equation).  The second
     16: !  equation  lambda = C * x enforces the surface flux auxiliary equation.  B and C have all
     17: !  positive entries, areas in C and fraction of area in B.
     18: !
     19: !!/*T
     20: !  Concepts: SNES^parallel Bratu example
     21: !  Concepts: MatNest
     22: !  Processors: n
     23: !T*/

     26: !
     27: !  --------------------------------------------------------------------------
     28: !
     29: !  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
     30: !  the partial differential equation
     31: !
     32: !          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
     33: !
     34: !  with boundary conditions
     35: !
     36: !           u = 0  for  x = 0, x = 1, y = 0, y = 1.
     37: !
     38: !  A finite difference approximation with the usual 5-point stencil
     39: !  is used to discretize the boundary value problem to obtain a nonlinear
     40: !  system of equations.
     41: !
     42: !  --------------------------------------------------------------------------
     43: !  The following define must be used before including any PETSc include files
     44: !  into a module or interface. This is because they can't handle declarations
     45: !  in them
     46: !
     47:       module petsc_kkt_solver
     48: #include <petsc/finclude/petscdm.h>
     49: #include <petsc/finclude/petscmat.h>
     50:       use petscdm
     51:       use petscmat
     52:       type petsc_kkt_solver_type
     53:         DM::da
     54: !     temp A block stuff
     55:         PetscInt mx,my
     56:         PetscMPIInt rank
     57:         PetscReal lambda
     58: !     Mats
     59:         Mat::Amat,AmatLin,Bmat,CMat,Dmat
     60:         IS::isPhi,isLambda
     61:       end type petsc_kkt_solver_type

     63:       end module petsc_kkt_solver

     65:       module petsc_kkt_solver_interfaces
     66:         use petsc_kkt_solver

     68:       Interface SNESSetApplicationContext
     69:         Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
     70: #include <petsc/finclude/petscsnes.h>
     71:         use petscsnes
     72:         use petsc_kkt_solver
     73:           SNES::    snesIn
     74:           type(petsc_kkt_solver_type) ctx
     75:           PetscErrorCode ierr
     76:         End Subroutine
     77:       End Interface SNESSetApplicationContext

     79:       Interface SNESGetApplicationContext
     80:         Subroutine SNESGetApplicationContext(snesIn,ctx,ierr)
     81: #include <petsc/finclude/petscsnes.h>
     82:         use petscsnes
     83:         use petsc_kkt_solver
     84:           SNES::     snesIn
     85:           type(petsc_kkt_solver_type), pointer :: ctx
     86:           PetscErrorCode ierr
     87:         End Subroutine
     88:       End Interface SNESGetApplicationContext
     89:       end module petsc_kkt_solver_interfaces

     91:       program main
     92: #include <petsc/finclude/petscdm.h>
     93: #include <petsc/finclude/petscsnes.h>
     94:       use petscdm
     95:       use petscdmda
     96:       use petscsnes
     97:       use petsc_kkt_solver
     98:       use petsc_kkt_solver_interfaces
     99:       implicit none
    100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101: !                   Variable declarations
    102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    103: !
    104: !  Variables:
    105: !     mysnes      - nonlinear solver
    106: !     x, r        - solution, residual vectors
    107: !     J           - Jacobian matrix
    108: !     its         - iterations for convergence
    109: !     Nx, Ny      - number of preocessors in x- and y- directions
    110: !
    111:       SNES::       mysnes
    112:       Vec::        x,r,x2,x1,x1loc,x2loc
    113:       Mat::       Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
    114: !      Mat::       tmat
    115:       DM::       daphi,dalam
    116:       IS::        isglobal(2)
    117:       PetscErrorCode   ierr
    118:       PetscInt         its,N1,N2,i,j,irow,row(1)
    119:       PetscInt         col(1),low,high,lamlow,lamhigh
    120:       PetscBool        flg
    121:       PetscInt         ione,nfour,itwo,nloc,nloclam
    122:       PetscReal lambda_max,lambda_min
    123:       type(petsc_kkt_solver_type)  solver
    124:       PetscScalar      bval(1),cval(1),one

    126: !  Note: Any user-defined Fortran routines (such as FormJacobian)
    127: !  MUST be declared as external.
    128:       external FormInitialGuess,FormJacobian,FormFunction

    130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131: !  Initialize program
    132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    133:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
    134:      if (ierr .ne. 0) then
    135:          print*,'PetscInitialize failed'
    136:          stop
    137:       endif
    138:       call MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr);CHKERRA(ierr)

    140: !  Initialize problem parameters
    141:       lambda_max  = 6.81_PETSC_REAL_KIND
    142:       lambda_min  = 0.0
    143:       solver%lambda = 6.0
    144:       ione = 1
    145:       nfour = 4
    146:       itwo = 2
    147:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr)
    148:       if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range'); endif

    150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    151: !  Create vector data structures; set function evaluation routine
    152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    154: !     just get size
    155:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione, &
    156:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,daphi,ierr);CHKERRA(ierr)
    157:       call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)
    158:       call DMSetUp(daphi,ierr);CHKERRA(ierr)
    159:       call DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                        &
    160:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
    161:       N1 = solver%my*solver%mx
    162:       N2 = solver%my
    163:       flg = .false.
    164:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr);CHKERRA(ierr)
    165:       if (flg) then
    166:          N2 = 0
    167:       endif

    169:       call DMDestroy(daphi,ierr);CHKERRA(ierr)

    171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172: !  Create matrix data structure; set Jacobian evaluation routine
    173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    174:       call DMShellCreate(PETSC_COMM_WORLD,daphi,ierr);CHKERRA(ierr)
    175:       call DMSetOptionsPrefix(daphi,'phi_',ierr);CHKERRA(ierr)
    176:       call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)

    178:       call VecCreate(PETSC_COMM_WORLD,x1,ierr);CHKERRA(ierr)
    179:       call VecSetSizes(x1,PETSC_DECIDE,N1,ierr);CHKERRA(ierr)
    180:       call VecSetFromOptions(x1,ierr);CHKERRA(ierr)

    182:       call VecGetOwnershipRange(x1,low,high,ierr);CHKERRA(ierr)
    183:       nloc = high - low

    185:       call MatCreate(PETSC_COMM_WORLD,Amat,ierr);CHKERRA(ierr)
    186:       call MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
    187:       call MatSetUp(Amat,ierr);CHKERRA(ierr)

    189:       call MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr);CHKERRA(ierr)
    190:       call MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
    191:       call MatSetUp(solver%AmatLin,ierr);CHKERRA(ierr)

    193:       call FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr);CHKERRA(ierr)
    194:       call MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    195:       call MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

    197:       call DMShellSetGlobalVector(daphi,x1,ierr);CHKERRA(ierr)
    198:       call DMShellSetMatrix(daphi,Amat,ierr);CHKERRA(ierr)

    200:       call VecCreate(PETSC_COMM_SELF,x1loc,ierr);CHKERRA(ierr)
    201:       call VecSetSizes(x1loc,nloc,nloc,ierr);CHKERRA(ierr)
    202:       call VecSetFromOptions(x1loc,ierr);CHKERRA(ierr)
    203:       call DMShellSetLocalVector(daphi,x1loc,ierr);CHKERRA(ierr)

    205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    206: !  Create B, C, & D matrices
    207: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    208:       call MatCreate(PETSC_COMM_WORLD,Cmat,ierr);CHKERRA(ierr)
    209:       call MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr);CHKERRA(ierr)
    210:       call MatSetUp(Cmat,ierr);CHKERRA(ierr)
    211: !      create data for C and B
    212:       call MatCreate(PETSC_COMM_WORLD,Bmat,ierr);CHKERRA(ierr)
    213:       call MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr);CHKERRA(ierr)
    214:       call MatSetUp(Bmat,ierr);CHKERRA(ierr)
    215: !     create data for D
    216:       call MatCreate(PETSC_COMM_WORLD,Dmat,ierr);CHKERRA(ierr)
    217:       call MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr);CHKERRA(ierr)
    218:       call MatSetUp(Dmat,ierr);CHKERRA(ierr)

    220:       call VecCreate(PETSC_COMM_WORLD,x2,ierr);CHKERRA(ierr)
    221:       call VecSetSizes(x2,PETSC_DECIDE,N2,ierr);CHKERRA(ierr)
    222:       call VecSetFromOptions(x2,ierr);CHKERRA(ierr)

    224:       call VecGetOwnershipRange(x2,lamlow,lamhigh,ierr);CHKERRA(ierr)
    225:       nloclam = lamhigh-lamlow

    227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    228: !  Set fake B and C
    229: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    230:       one    = 1.0
    231:       if (N2 .gt. 0) then
    232:          bval(1) = -one/(solver%mx-2)
    233: !     cval = -one/(solver%my*solver%mx)
    234:          cval(1) = -one
    235:          do 20 irow=low,high-1
    236:             j = irow/solver%mx   ! row in domain
    237:             i = mod(irow,solver%mx)
    238:             row(1) = irow
    239:             col(1) = j
    240:             if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
    241:                !     no op
    242:             else
    243:                call MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr);CHKERRA(ierr)
    244:             endif
    245:             row(1) = j
    246:             call MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
    247:  20   continue
    248:       endif
    249:       call MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    250:       call MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    251:       call MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    252:       call MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

    254: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    255: !  Set D (indentity)
    256: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    257:       do 30 j=lamlow,lamhigh-1
    258:          row(1) = j
    259:          cval(1) = one
    260:          call MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
    261:  30   continue
    262:       call MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    263:       call MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

    265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    266: !  DM for lambda (dalam) : temp driver for A block, setup A block solver data
    267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    268:       call DMShellCreate(PETSC_COMM_WORLD,dalam,ierr);CHKERRA(ierr)
    269:       call DMShellSetGlobalVector(dalam,x2,ierr);CHKERRA(ierr)
    270:       call DMShellSetMatrix(dalam,Dmat,ierr);CHKERRA(ierr)

    272:       call VecCreate(PETSC_COMM_SELF,x2loc,ierr);CHKERRA(ierr)
    273:       call VecSetSizes(x2loc,nloclam,nloclam,ierr);CHKERRA(ierr)
    274:       call VecSetFromOptions(x2loc,ierr);CHKERRA(ierr)
    275:       call DMShellSetLocalVector(dalam,x2loc,ierr);CHKERRA(ierr)

    277:       call DMSetOptionsPrefix(dalam,'lambda_',ierr);CHKERRA(ierr)
    278:       call DMSetFromOptions(dalam,ierr);CHKERRA(ierr)
    279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    280: !  Create field split DA
    281: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    282:       call DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr);CHKERRA(ierr)
    283:       call DMCompositeAddDM(solver%da,daphi,ierr);CHKERRA(ierr)
    284:       call DMCompositeAddDM(solver%da,dalam,ierr);CHKERRA(ierr)
    285:       call DMSetFromOptions(solver%da,ierr);CHKERRA(ierr)
    286:       call DMSetUp(solver%da,ierr);CHKERRA(ierr)
    287:       call DMCompositeGetGlobalISs(solver%da,isglobal,ierr);CHKERRA(ierr)
    288:       solver%isPhi = isglobal(1)
    289:       solver%isLambda = isglobal(2)

    291: !     cache matrices
    292:       solver%Amat = Amat
    293:       solver%Bmat = Bmat
    294:       solver%Cmat = Cmat
    295:       solver%Dmat = Dmat

    297:       matArray(1) = Amat
    298:       matArray(2) = Bmat
    299:       matArray(3) = Cmat
    300:       matArray(4) = Dmat

    302:       call MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr);CHKERRA(ierr)
    303:       call MatSetFromOptions(KKTmat,ierr);CHKERRA(ierr)

    305: !  Extract global and local vectors from DMDA; then duplicate for remaining
    306: !     vectors that are the same types
    307:       call MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
    308:       call VecDuplicate(x,r,ierr);CHKERRA(ierr)

    310:       call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr)

    312:       call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)

    314:       call SNESSetApplicationContext(mysnes,solver,ierr);CHKERRA(ierr)

    316:       call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)

    318: !  Set function evaluation routine and vector
    319:       call SNESSetFunction(mysnes,r,FormFunction,solver,ierr);CHKERRA(ierr)

    321:       call SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr);CHKERRA(ierr)

    323: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    324: !  Customize nonlinear solver; set runtime options
    325: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    326: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
    327:       call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr)

    329: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    330: !  Evaluate initial guess; then solve nonlinear system.
    331: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    332: !  Note: The user should initialize the vector, x, with the initial guess
    333: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
    334: !  to employ an initial guess of zero, the user should explicitly set
    335: !  this vector to zero by calling VecSet().

    337:       call FormInitialGuess(mysnes,x,ierr);CHKERRA(ierr)
    338:       call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr)
    339:       call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr)
    340:       if (solver%rank .eq. 0) then
    341:          write(6,100) its
    342:       endif
    343:   100 format('Number of SNES iterations = ',i5)

    345: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    346: !  Free work space.  All PETSc objects should be destroyed when they
    347: !  are no longer needed.
    348: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    349:       call MatDestroy(KKTmat,ierr);CHKERRA(ierr)
    350:       call MatDestroy(Amat,ierr);CHKERRA(ierr)
    351:       call MatDestroy(Dmat,ierr);CHKERRA(ierr)
    352:       call MatDestroy(Bmat,ierr);CHKERRA(ierr)
    353:       call MatDestroy(Cmat,ierr);CHKERRA(ierr)
    354:       call MatDestroy(solver%AmatLin,ierr);CHKERRA(ierr)
    355:       call ISDestroy(solver%isPhi,ierr);CHKERRA(ierr)
    356:       call ISDestroy(solver%isLambda,ierr);CHKERRA(ierr)
    357:       call VecDestroy(x,ierr);CHKERRA(ierr)
    358:       call VecDestroy(x2,ierr);CHKERRA(ierr)
    359:       call VecDestroy(x1,ierr);CHKERRA(ierr)
    360:       call VecDestroy(x1loc,ierr);CHKERRA(ierr)
    361:       call VecDestroy(x2loc,ierr);CHKERRA(ierr)
    362:       call VecDestroy(r,ierr);CHKERRA(ierr)
    363:       call SNESDestroy(mysnes,ierr);CHKERRA(ierr)
    364:       call DMDestroy(solver%da,ierr);CHKERRA(ierr)
    365:       call DMDestroy(daphi,ierr);CHKERRA(ierr)
    366:       call DMDestroy(dalam,ierr);CHKERRA(ierr)

    368:       call PetscFinalize(ierr)
    369:       end

    371: ! ---------------------------------------------------------------------
    372: !
    373: !  FormInitialGuess - Forms initial approximation.
    374: !
    375: !  Input Parameters:
    376: !  X - vector
    377: !
    378: !  Output Parameter:
    379: !  X - vector
    380: !
    381: !  Notes:
    382: !  This routine serves as a wrapper for the lower-level routine
    383: !  "InitialGuessLocal", where the actual computations are
    384: !  done using the standard Fortran style of treating the local
    385: !  vector data as a multidimensional array over the local mesh.
    386: !  This routine merely handles ghost point scatters and accesses
    387: !  the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
    388: !
    389:       subroutine FormInitialGuess(mysnes,Xnest,ierr)
    390: #include <petsc/finclude/petscsnes.h>
    391:       use petscsnes
    392:       use petsc_kkt_solver
    393:       use petsc_kkt_solver_interfaces
    394:       implicit none
    395: !  Input/output variables:
    396:       SNES::     mysnes
    397:       Vec::      Xnest
    398:       PetscErrorCode ierr

    400: !  Declarations for use with local arrays:
    401:       type(petsc_kkt_solver_type), pointer:: solver
    402:       Vec::      Xsub(2)
    403:       PetscInt::  izero,ione,itwo

    405:       izero = 0
    406:       ione = 1
    407:       itwo = 2
    408:       0
    409:       call SNESGetApplicationContext(mysnes,solver,ierr);CHKERRQ(ierr)
    410:       call DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

    412:       call InitialGuessLocal(solver,Xsub(1),ierr);CHKERRQ(ierr)
    413:       call VecAssemblyBegin(Xsub(1),ierr);CHKERRQ(ierr)
    414:       call VecAssemblyEnd(Xsub(1),ierr);CHKERRQ(ierr)

    416: !     zero out lambda
    417:       call VecZeroEntries(Xsub(2),ierr);CHKERRQ(ierr)
    418:       call DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

    420:       return
    421:       end subroutine FormInitialGuess

    423: ! ---------------------------------------------------------------------
    424: !
    425: !  InitialGuessLocal - Computes initial approximation, called by
    426: !  the higher level routine FormInitialGuess().
    427: !
    428: !  Input Parameter:
    429: !  X1 - local vector data
    430: !
    431: !  Output Parameters:
    432: !  x - local vector data
    433: !  ierr - error code
    434: !
    435: !  Notes:
    436: !  This routine uses standard Fortran-style computations over a 2-dim array.
    437: !
    438:       subroutine InitialGuessLocal(solver,X1,ierr)
    439: #include <petsc/finclude/petscsys.h>
    440:       use petscsys
    441:       use petsc_kkt_solver
    442:       implicit none
    443: !  Input/output variables:
    444:       type (petsc_kkt_solver_type)         solver
    445:       Vec::      X1
    446:       PetscErrorCode ierr

    448: !  Local variables:
    449:       PetscInt      row,i,j,ione,low,high
    450:       PetscReal   temp1,temp,hx,hy,v
    451:       PetscReal   one

    453: !  Set parameters
    454:       ione = 1
    455:       0
    456:       one    = 1.0
    457:       hx     = one/(solver%mx-1)
    458:       hy     = one/(solver%my-1)
    459:       temp1  = solver%lambda/(solver%lambda + one) + one

    461:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)

    463:       do 20 row=low,high-1
    464:          j = row/solver%mx
    465:          i = mod(row,solver%mx)
    466:          temp = min(j,solver%my-j+1)*hy
    467:          if (i .eq. 0 .or. j .eq. 0  .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
    468:             v = 0.0
    469:          else
    470:             v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
    471:          endif
    472:          call VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
    473:  20   continue

    475:       return
    476:       end subroutine InitialGuessLocal

    479: ! ---------------------------------------------------------------------
    480: !
    481: !  FormJacobian - Evaluates Jacobian matrix.
    482: !
    483: !  Input Parameters:
    484: !  dummy     - the SNES context
    485: !  x         - input vector
    486: !  solver    - solver data
    487: !
    488: !  Output Parameters:
    489: !  jac      - Jacobian matrix
    490: !  jac_prec - optionally different preconditioning matrix (not used here)
    491: !  flag     - flag indicating matrix structure
    492: !
    493:       subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
    494: #include <petsc/finclude/petscsnes.h>
    495:       use petscsnes
    496:       use petsc_kkt_solver
    497:       implicit none
    498: !  Input/output variables:
    499:       SNES::     dummy
    500:       Vec::      X
    501:      Mat::     jac,jac_prec
    502:       type(petsc_kkt_solver_type)  solver
    503:       PetscErrorCode ierr

    505: !  Declarations for use with local arrays:
    506:       Vec::      Xsub(1)
    507:      Mat::     Amat
    508:       PetscInt       izero,ione

    510:       izero = 0
    511:       ione = 1

    513:       call DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

    515: !     Compute entries for the locally owned part of the Jacobian preconditioner.
    516:       call MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr);CHKERRQ(ierr)

    518:       call FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr);CHKERRQ(ierr)
    519:       call MatDestroy(Amat,ierr);CHKERRQ(ierr) ! discard our reference
    520:       call DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

    522:       ! the rest of the matrix is not touched
    523:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    524:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    525:       if (jac .ne. jac_prec) then
    526:          call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    527:          call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    528:       end if

    530: !     Tell the matrix we will never add a new nonzero location to the
    531: !     matrix. If we do it will generate an error.
    532:       call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr);CHKERRQ(ierr)

    534:       return
    535:       end subroutine FormJacobian

    537: ! ---------------------------------------------------------------------
    538: !
    539: !  FormJacobianLocal - Computes Jacobian preconditioner matrix,
    540: !  called by the higher level routine FormJacobian().
    541: !
    542: !  Input Parameters:
    543: !  x        - local vector data
    544: !
    545: !  Output Parameters:
    546: !  jac - Jacobian preconditioner matrix
    547: !  ierr     - error code
    548: !
    549: !  Notes:
    550: !  This routine uses standard Fortran-style computations over a 2-dim array.
    551: !
    552:       subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
    553: #include <petsc/finclude/petscmat.h>
    554:       use petscmat
    555:       use petsc_kkt_solver
    556:       implicit none
    557: !  Input/output variables:
    558:       type (petsc_kkt_solver_type) solver
    559:       Vec::      X1
    560:      Mat::     jac
    561:       logical        add_nl_term
    562:       PetscErrorCode ierr

    564: !  Local variables:
    565:       PetscInt    irow,row(1),col(5),i,j
    566:       PetscInt    ione,ifive,low,high,ii
    567:       PetscScalar two,one,hx,hy,hy2inv
    568:       PetscScalar hx2inv,sc,v(5)
    569:       PetscScalar,pointer :: lx_v(:)

    571: !  Set parameters
    572:       ione   = 1
    573:       ifive  = 5
    574:       one    = 1.0
    575:       two    = 2.0
    576:       hx     = one/(solver%mx-1)
    577:       hy     = one/(solver%my-1)
    578:       sc     = solver%lambda
    579:       hx2inv = one/(hx*hx)
    580:       hy2inv = one/(hy*hy)

    582:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
    583:       call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

    585:       ii = 0
    586:       do 20 irow=low,high-1
    587:          j = irow/solver%mx
    588:          i = mod(irow,solver%mx)
    589:          ii = ii + 1            ! one based local index
    590: !     boundary points
    591:          if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
    592:             col(1) = irow
    593:             row(1) = irow
    594:             v(1)   = one
    595:             call MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
    596: !     interior grid points
    597:          else
    598:             v(1) = -hy2inv
    599:             if (j-1==0) v(1) = 0.0
    600:             v(2) = -hx2inv
    601:             if (i-1==0) v(2) = 0.0
    602:             v(3) = two*(hx2inv + hy2inv)
    603:             if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
    604:             v(4) = -hx2inv
    605:             if (i+1==solver%mx-1) v(4) = 0.0
    606:             v(5) = -hy2inv
    607:             if (j+1==solver%my-1) v(5) = 0.0
    608:             col(1) = irow - solver%mx
    609:             col(2) = irow - 1
    610:             col(3) = irow
    611:             col(4) = irow + 1
    612:             col(5) = irow + solver%mx
    613:             row(1) = irow
    614:             call MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
    615:          endif
    616:  20   continue

    618:       call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

    620:       return
    621:       end subroutine FormJacobianLocal

    624: ! ---------------------------------------------------------------------
    625: !
    626: !  FormFunction - Evaluates nonlinear function, F(x).
    627: !
    628: !  Input Parameters:
    629: !  snes - the SNES context
    630: !  X - input vector
    631: !  dummy - optional user-defined context, as set by SNESSetFunction()
    632: !          (not used here)
    633: !
    634: !  Output Parameter:
    635: !  F - function vector
    636: !
    637:       subroutine FormFunction(snesIn,X,F,solver,ierr)
    638: #include <petsc/finclude/petscsnes.h>
    639:       use petscsnes
    640:       use petsc_kkt_solver
    641:       implicit none
    642: !  Input/output variables:
    643:       SNES::     snesIn
    644:      Vec::      X,F
    645:       PetscErrorCode ierr
    646:       type (petsc_kkt_solver_type) solver

    648: !  Declarations for use with local arrays:
    649:      Vec::              Xsub(2),Fsub(2)
    650:       PetscInt               izero,ione,itwo

    652: !  Scatter ghost points to local vector, using the 2-step process
    653: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
    654: !  By placing code between these two statements, computations can
    655: !  be done while messages are in transition.

    657:       izero = 0
    658:       ione = 1
    659:       itwo = 2
    660:       call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
    661:       call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)

    663:       call FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr);CHKERRQ(ierr)
    664:       call MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)

    666: !     do rest of operator (linear)
    667:       call MatMult(    solver%Cmat, Xsub(1),      Fsub(2), ierr);CHKERRQ(ierr)
    668:       call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
    669:       call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr)

    671:       call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
    672:       call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
    673:       return
    674:       end subroutine formfunction

    677: ! ---------------------------------------------------------------------
    678: !
    679: !  FormFunctionNLTerm - Computes nonlinear function, called by
    680: !  the higher level routine FormFunction().
    681: !
    682: !  Input Parameter:
    683: !  x - local vector data
    684: !
    685: !  Output Parameters:
    686: !  f - local vector data, f(x)
    687: !  ierr - error code
    688: !
    689: !  Notes:
    690: !  This routine uses standard Fortran-style computations over a 2-dim array.
    691: !
    692:       subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
    693: #include <petsc/finclude/petscvec.h>
    694:       use petscvec
    695:       use petsc_kkt_solver
    696:       implicit none
    697: !  Input/output variables:
    698:       type (petsc_kkt_solver_type) solver
    699:      Vec::      X1,F1
    700:       PetscErrorCode ierr
    701: !  Local variables:
    702:       PetscScalar one,sc
    703:       PetscScalar u,v(1)
    704:       PetscInt  i,j,low,high,ii,ione,irow,row(1)
    705:       PetscScalar,pointer :: lx_v(:)

    707:       one    = 1.0
    708:       sc     = solver%lambda
    709:       ione   = 1

    711:       call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
    712:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)

    714: !     Compute function over the locally owned part of the grid
    715:       ii = 0
    716:       do 20 irow=low,high-1
    717:          j = irow/solver%mx
    718:          i = mod(irow,solver%mx)
    719:          ii = ii + 1            ! one based local index
    720:          row(1) = irow
    721:          if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
    722:             v(1) = 0.0
    723:          else
    724:             u = lx_v(ii)
    725:             v(1) = -sc*exp(u)
    726:          endif
    727:          call VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
    728:  20   continue

    730:       call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

    732:       call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr)
    733:       call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr)

    735:       0
    736:       return
    737:       end subroutine FormFunctionNLTerm

    739: !/*TEST
    740: !
    741: !   build:
    742: !      requires: !single !complex
    743: !
    744: !   test:
    745: !      nsize: 4
    746: !      args: -par 5.0 -da_grid_x 10 -da_grid_y 10 -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason -ksp_type fgmres -ksp_norm_type unpreconditioned -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -ksp_monitor_short -fieldsplit_lambda_ksp_type preonly -fieldsplit_lambda_pc_type jacobi -fieldsplit_phi_pc_type gamg -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0. -fieldsplit_phi_pc_gamg_esteig_ksp_type cg
    747: !
    748: !TEST*/