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: !

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

 60:   interface
 61:     subroutine SNESSetApplicationContext(snesIn, ctx, ierr)
 62:       use petscsnes
 63:       import ex73f90tmodule_type
 64:       SNES :: snesIn
 65:       type(ex73f90tmodule_type) ctx
 66:       PetscErrorCode ierr
 67:     end subroutine
 68:     subroutine SNESGetApplicationContext(snesIn, ctx, ierr)
 69:       use petscsnes
 70:       import ex73f90tmodule_type
 71:       SNES :: snesIn
 72:       type(ex73f90tmodule_type), pointer :: ctx
 73:       PetscErrorCode ierr
 74:     end subroutine
 75:   end interface

 77: contains
 78:   subroutine MyObjective(snes, x, result, ctx, ierr)
 79:     PetscInt ctx
 80:     Vec x, f
 81:     SNES snes
 82:     PetscErrorCode ierr
 83:     PetscScalar result
 84:     PetscReal fnorm

 86:     PetscCall(VecDuplicate(x, f, ierr))
 87:     PetscCall(SNESComputeFunction(snes, x, f, ierr))
 88:     PetscCall(VecNorm(f, NORM_2, fnorm, ierr))
 89:     result = .5*fnorm*fnorm
 90:     PetscCall(VecDestroy(f, ierr))
 91:   end subroutine MyObjective

 93: ! ---------------------------------------------------------------------
 94: !
 95: !  FormInitialGuess - Forms initial approximation.
 96: !
 97: !  Input Parameters:
 98: !  X - vector
 99: !
100: !  Output Parameter:
101: !  X - vector
102: !
103: !  Notes:
104: !  This routine serves as a wrapper for the lower-level routine
105: !  "InitialGuessLocal", where the actual computations are
106: !  done using the standard Fortran style of treating the local
107: !  vector data as a multidimensional array over the local mesh.
108: !  This routine merely handles ghost point scatters and accesses
109: !  the local vector data via VecGetArray() and VecRestoreArray().
110: !
111:   subroutine FormInitialGuess(mysnes, Xnest, ierr)
112: !  Input/output variables:
113:     SNES::     mysnes
114:     Vec::      Xnest
115:     PetscErrorCode ierr

117: !  Declarations for use with local arrays:
118:     type(ex73f90tmodule_type), pointer:: solver
119:     Vec::      Xsub(2)
120:     PetscInt::  izero, ione, itwo

122:     izero = 0
123:     ione = 1
124:     itwo = 2
125:     ierr = 0
126:     PetscCall(SNESGetApplicationContext(mysnes, solver, ierr))
127:     PetscCall(DMCompositeGetAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

129:     PetscCall(InitialGuessLocal(solver, Xsub(1), ierr))
130:     PetscCall(VecAssemblyBegin(Xsub(1), ierr))
131:     PetscCall(VecAssemblyEnd(Xsub(1), ierr))

133: !     zero out lambda
134:     PetscCall(VecZeroEntries(Xsub(2), ierr))
135:     PetscCall(DMCompositeRestoreAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

137:   end subroutine FormInitialGuess

139: ! ---------------------------------------------------------------------
140: !
141: !  InitialGuessLocal - Computes initial approximation, called by
142: !  the higher level routine FormInitialGuess().
143: !
144: !  Input Parameter:
145: !  X1 - local vector data
146: !
147: !  Output Parameters:
148: !  x - local vector data
149: !  ierr - error code
150: !
151: !  Notes:
152: !  This routine uses standard Fortran-style computations over a 2-dim array.
153: !
154:   subroutine InitialGuessLocal(solver, X1, ierr)
155: !  Input/output variables:
156:     type(ex73f90tmodule_type) solver
157:     Vec::      X1
158:     PetscErrorCode ierr

160: !  Local variables:
161:     PetscInt row, i, j, ione, low, high
162:     PetscReal temp1, temp, hx, hy, v
163:     PetscReal one

165: !  Set parameters
166:     ione = 1
167:     ierr = 0
168:     one = 1.0
169:     hx = one/(solver%mx - 1)
170:     hy = one/(solver%my - 1)
171:     temp1 = solver%lambda/(solver%lambda + one) + one

173:     PetscCall(VecGetOwnershipRange(X1, low, high, ierr))

175:     do row = low, high - 1
176:       j = row/solver%mx
177:       i = mod(row, solver%mx)
178:       temp = min(j, solver%my - j + 1)*hy
179:       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
180:         v = 0.0
181:       else
182:         v = temp1*sqrt(min(min(i, solver%mx - i + 1)*hx, temp))
183:       end if
184:       PetscCall(VecSetValues(X1, ione, [row], [v], INSERT_VALUES, ierr))
185:     end do

187:   end subroutine InitialGuessLocal

189: ! ---------------------------------------------------------------------
190: !
191: !  FormJacobian - Evaluates Jacobian matrix.
192: !
193: !  Input Parameters:
194: !  dummy     - the SNES context
195: !  x         - input vector
196: !  solver    - solver data
197: !
198: !  Output Parameters:
199: !  jac      - Jacobian matrix
200: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
201: !
202:   subroutine FormJacobian(dummy, X, jac, jac_prec, solver, ierr)
203: !  Input/output variables:
204:     SNES::     dummy
205:     Vec::      X
206:     Mat::     jac, jac_prec
207:     type(ex73f90tmodule_type) solver
208:     PetscErrorCode ierr

210: !  Declarations for use with local arrays:
211:     Vec::      Xsub(1)
212:     Mat::     Amat
213:     PetscInt ione

215:     ione = 1

217:     PetscCall(DMCompositeGetAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

219: !     Compute entries for the locally owned part of the Jacobian preconditioner.
220:     PetscCall(MatCreateSubMatrix(jac_prec, solver%isPhi, solver%isPhi, MAT_INITIAL_MATRIX, Amat, ierr))

222:     PetscCall(FormJacobianLocal(Xsub(1), Amat, solver, .true., ierr))
223:     PetscCall(MatDestroy(Amat, ierr)) ! discard our reference
224:     PetscCall(DMCompositeRestoreAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

226:     ! the rest of the matrix is not touched
227:     PetscCall(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
228:     PetscCall(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
229:     if (jac /= jac_prec) then
230:       PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
231:       PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
232:     end if

234: !     Tell the matrix we will never add a new nonzero location to the
235: !     matrix. If we do it will generate an error.
236:     PetscCall(MatSetOption(jac_prec, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))

238:   end subroutine FormJacobian

240: ! ---------------------------------------------------------------------
241: !
242: !  FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
243: !  called by the higher level routine FormJacobian().
244: !
245: !  Input Parameters:
246: !  x        - local vector data
247: !
248: !  Output Parameters:
249: !  jac - Jacobian matrix used to compute the preconditioner
250: !  ierr     - error code
251: !
252: !  Notes:
253: !  This routine uses standard Fortran-style computations over a 2-dim array.
254: !
255:   subroutine FormJacobianLocal(X1, jac, solver, add_nl_term, ierr)
256: !  Input/output variables:
257:     type(ex73f90tmodule_type) solver
258:     Vec::      X1
259:     Mat::     jac
260:     logical add_nl_term
261:     PetscErrorCode ierr

263: !  Local variables:
264:     PetscInt irow, row(1), col(5), i, j
265:     PetscInt ione, ifive, low, high, ii
266:     PetscScalar two, one, hx, hy, hy2inv
267:     PetscScalar hx2inv, sc, v(5)
268:     PetscScalar, pointer :: lx_v(:)

270: !  Set parameters
271:     ione = 1
272:     ifive = 5
273:     one = 1.0
274:     two = 2.0
275:     hx = one/(solver%mx - 1)
276:     hy = one/(solver%my - 1)
277:     sc = solver%lambda
278:     hx2inv = one/(hx*hx)
279:     hy2inv = one/(hy*hy)

281:     PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
282:     PetscCall(VecGetArrayRead(X1, lx_v, ierr))

284:     ii = 0
285:     do irow = low, high - 1
286:       j = irow/solver%mx
287:       i = mod(irow, solver%mx)
288:       ii = ii + 1            ! one based local index
289: !     boundary points
290:       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
291:         col(1) = irow
292:         row(1) = irow
293:         v(1) = one
294:         PetscCall(MatSetValues(jac, ione, row, ione, col, v, INSERT_VALUES, ierr))
295: !     interior grid points
296:       else
297:         v(1) = -hy2inv
298:         if (j - 1 == 0) v(1) = 0.0
299:         v(2) = -hx2inv
300:         if (i - 1 == 0) v(2) = 0.0
301:         v(3) = two*(hx2inv + hy2inv)
302:         if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
303:         v(4) = -hx2inv
304:         if (i + 1 == solver%mx - 1) v(4) = 0.0
305:         v(5) = -hy2inv
306:         if (j + 1 == solver%my - 1) v(5) = 0.0
307:         col(1) = irow - solver%mx
308:         col(2) = irow - 1
309:         col(3) = irow
310:         col(4) = irow + 1
311:         col(5) = irow + solver%mx
312:         row(1) = irow
313:         PetscCall(MatSetValues(jac, ione, row, ifive, col, v, INSERT_VALUES, ierr))
314:       end if
315:     end do

317:     PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))

319:   end subroutine FormJacobianLocal

321: ! ---------------------------------------------------------------------
322: !
323: !  FormFunction - Evaluates nonlinear function, F(x).
324: !
325: !  Input Parameters:
326: !  snes - the SNES context
327: !  X - input vector
328: !  dummy - optional user-defined context, as set by SNESSetFunction()
329: !          (not used here)
330: !
331: !  Output Parameter:
332: !  F - function vector
333: !
334:   subroutine FormFunction(snesIn, X, F, solver, ierr)
335: !  Input/output variables:
336:     SNES::     snesIn
337:     Vec::      X, F
338:     PetscErrorCode ierr
339:     type(ex73f90tmodule_type) solver

341: !  Declarations for use with local arrays:
342:     Vec::              Xsub(2), Fsub(2)
343:     PetscInt itwo

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

350:     itwo = 2
351:     PetscCall(DMCompositeGetAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
352:     PetscCall(DMCompositeGetAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))

354:     PetscCall(FormFunctionNLTerm(Xsub(1), Fsub(1), solver, ierr))
355:     PetscCall(MatMultAdd(solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))

357: !     do rest of operator (linear)
358:     PetscCall(MatMult(solver%Cmat, Xsub(1), Fsub(2), ierr))
359:     PetscCall(MatMultAdd(solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
360:     PetscCall(MatMultAdd(solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))

362:     PetscCall(DMCompositeRestoreAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
363:     PetscCall(DMCompositeRestoreAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
364:   end subroutine formfunction

366: ! ---------------------------------------------------------------------
367: !
368: !  FormFunctionNLTerm - Computes nonlinear function, called by
369: !  the higher level routine FormFunction().
370: !
371: !  Input Parameter:
372: !  x - local vector data
373: !
374: !  Output Parameters:
375: !  f - local vector data, f(x)
376: !  ierr - error code
377: !
378: !  Notes:
379: !  This routine uses standard Fortran-style computations over a 2-dim array.
380: !
381:   subroutine FormFunctionNLTerm(X1, F1, solver, ierr)
382: !  Input/output variables:
383:     type(ex73f90tmodule_type) solver
384:     Vec::      X1, F1
385:     PetscErrorCode ierr
386: !  Local variables:
387:     PetscScalar sc
388:     PetscScalar u, v(1)
389:     PetscInt i, j, low, high, ii, ione, irow, row(1)
390:     PetscScalar, pointer :: lx_v(:)

392:     sc = solver%lambda
393:     ione = 1

395:     PetscCall(VecGetArrayRead(X1, lx_v, ierr))
396:     PetscCall(VecGetOwnershipRange(X1, low, high, ierr))

398: !     Compute function over the locally owned part of the grid
399:     ii = 0
400:     do irow = low, high - 1
401:       j = irow/solver%mx
402:       i = mod(irow, solver%mx)
403:       ii = ii + 1            ! one based local index
404:       row(1) = irow
405:       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
406:         v(1) = 0.0
407:       else
408:         u = lx_v(ii)
409:         v(1) = -sc*exp(u)
410:       end if
411:       PetscCall(VecSetValues(F1, ione, row, v, INSERT_VALUES, ierr))
412:     end do

414:     PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))

416:     PetscCall(VecAssemblyBegin(F1, ierr))
417:     PetscCall(VecAssemblyEnd(F1, ierr))

419:     ierr = 0
420:   end subroutine FormFunctionNLTerm

422: end module ex73f90tmodule

424: program main
425:   use petscsnes
426:   use ex73f90tmodule
427:   implicit none
428: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429: !                   Variable declarations
430: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431: !
432: !  Variables:
433: !     mysnes      - nonlinear solver
434: !     x, r        - solution, residual vectors
435: !     J           - Jacobian matrix
436: !     its         - iterations for convergence
437: !     Nx, Ny      - number of preocessors in x- and y- directions
438: !
439:   SNES::       mysnes
440:   Vec::        x, r, x2, x1, x1loc, x2loc
441:   Mat::       Amat, Bmat, Cmat, Dmat, KKTMat, matArray(4)
442: !      Mat::       tmat
443:   DM::       daphi, dalam
444:   IS, pointer ::        isglobal(:)
445:   PetscErrorCode ierr
446:   PetscInt its, N1, N2, i, j, irow, row(1)
447:   PetscInt col(1), low, high, lamlow, lamhigh
448:   PetscBool flg
449:   PetscInt ione, nfour, itwo, nloc, nloclam
450:   PetscReal lambda_max, lambda_min
451:   type(ex73f90tmodule_type) solver
452:   PetscScalar bval(1), cval(1), one
453:   PetscBool useobjective

455: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456: !  Initialize program
457: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
458:   PetscCallA(PetscInitialize(ierr))
459:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, solver%rank, ierr))

461: !  Initialize problem parameters
462:   lambda_max = 6.81_PETSC_REAL_KIND
463:   lambda_min = 0.0
464:   solver%lambda = 6.0
465:   ione = 1
466:   nfour = 4
467:   itwo = 2
468:   useobjective = PETSC_FALSE
469:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', solver%lambda, flg, ierr))
470:   PetscCheckA(solver%lambda <= lambda_max .and. solver%lambda >= lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')
471:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-objective', useobjective, flg, ierr))

473: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
474: !  Create vector data structures; set function evaluation routine
475: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

477: !     just get size
478:   PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nfour, nfour, PETSC_DECIDE, PETSC_DECIDE, ione, ione, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, daphi, ierr))
479:   PetscCallA(DMSetFromOptions(daphi, ierr))
480:   PetscCallA(DMSetUp(daphi, ierr))
481:   PetscCallA(DMDAGetInfo(daphi, PETSC_NULL_INTEGER, solver%mx, solver%my, PETSC_NULL_INTEGER, 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))
482:   N1 = solver%my*solver%mx
483:   N2 = solver%my
484:   flg = .false.
485:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_constraints', flg, flg, ierr))
486:   if (flg) then
487:     N2 = 0
488:   end if

490:   PetscCallA(DMDestroy(daphi, ierr))

492: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
493: !  Create matrix data structure; set Jacobian evaluation routine
494: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
495:   PetscCallA(DMShellCreate(PETSC_COMM_WORLD, daphi, ierr))
496:   PetscCallA(DMSetOptionsPrefix(daphi, 'phi_', ierr))
497:   PetscCallA(DMSetFromOptions(daphi, ierr))

499:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x1, ierr))
500:   PetscCallA(VecSetSizes(x1, PETSC_DECIDE, N1, ierr))
501:   PetscCallA(VecSetFromOptions(x1, ierr))

503:   PetscCallA(VecGetOwnershipRange(x1, low, high, ierr))
504:   nloc = high - low

506:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr))
507:   PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
508:   PetscCallA(MatSetUp(Amat, ierr))

510:   PetscCallA(MatCreate(PETSC_COMM_WORLD, solver%AmatLin, ierr))
511:   PetscCallA(MatSetSizes(solver%AmatLin, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
512:   PetscCallA(MatSetUp(solver%AmatLin, ierr))

514:   PetscCallA(FormJacobianLocal(x1, solver%AmatLin, solver, .false., ierr))
515:   PetscCallA(MatAssemblyBegin(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
516:   PetscCallA(MatAssemblyEnd(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))

518:   PetscCallA(DMShellSetGlobalVector(daphi, x1, ierr))
519:   PetscCallA(DMShellSetMatrix(daphi, Amat, ierr))

521:   PetscCallA(VecCreate(PETSC_COMM_SELF, x1loc, ierr))
522:   PetscCallA(VecSetSizes(x1loc, nloc, nloc, ierr))
523:   PetscCallA(VecSetFromOptions(x1loc, ierr))
524:   PetscCallA(DMShellSetLocalVector(daphi, x1loc, ierr))

526: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
527: !  Create B, C, & D matrices
528: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
529:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Cmat, ierr))
530:   PetscCallA(MatSetSizes(Cmat, PETSC_DECIDE, PETSC_DECIDE, N2, N1, ierr))
531:   PetscCallA(MatSetUp(Cmat, ierr))
532: !      create data for C and B
533:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Bmat, ierr))
534:   PetscCallA(MatSetSizes(Bmat, PETSC_DECIDE, PETSC_DECIDE, N1, N2, ierr))
535:   PetscCallA(MatSetUp(Bmat, ierr))
536: !     create data for D
537:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Dmat, ierr))
538:   PetscCallA(MatSetSizes(Dmat, PETSC_DECIDE, PETSC_DECIDE, N2, N2, ierr))
539:   PetscCallA(MatSetUp(Dmat, ierr))

541:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x2, ierr))
542:   PetscCallA(VecSetSizes(x2, PETSC_DECIDE, N2, ierr))
543:   PetscCallA(VecSetFromOptions(x2, ierr))

545:   PetscCallA(VecGetOwnershipRange(x2, lamlow, lamhigh, ierr))
546:   nloclam = lamhigh - lamlow

548: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
549: !  Set fake B and C
550: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
551:   one = 1.0
552:   if (N2 > 0) then
553:     bval(1) = -one/(solver%mx - 2)
554: !     cval = -one/(solver%my*solver%mx)
555:     cval(1) = -one
556:     do irow = low, high - 1
557:       j = irow/solver%mx   ! row in domain
558:       i = mod(irow, solver%mx)
559:       row(1) = irow
560:       col(1) = j
561:       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
562:         !     no op
563:       else
564:         PetscCallA(MatSetValues(Bmat, ione, row, ione, col, bval, INSERT_VALUES, ierr))
565:       end if
566:       row(1) = j
567:       PetscCallA(MatSetValues(Cmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
568:     end do
569:   end if
570:   PetscCallA(MatAssemblyBegin(Bmat, MAT_FINAL_ASSEMBLY, ierr))
571:   PetscCallA(MatAssemblyEnd(Bmat, MAT_FINAL_ASSEMBLY, ierr))
572:   PetscCallA(MatAssemblyBegin(Cmat, MAT_FINAL_ASSEMBLY, ierr))
573:   PetscCallA(MatAssemblyEnd(Cmat, MAT_FINAL_ASSEMBLY, ierr))

575: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
576: !  Set D (identity)
577: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578:   do j = lamlow, lamhigh - 1
579:     row(1) = j
580:     cval(1) = one
581:     PetscCallA(MatSetValues(Dmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
582:   end do
583:   PetscCallA(MatAssemblyBegin(Dmat, MAT_FINAL_ASSEMBLY, ierr))
584:   PetscCallA(MatAssemblyEnd(Dmat, MAT_FINAL_ASSEMBLY, ierr))

586: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
587: !  DM for lambda (dalam) : temp driver for A block, setup A block solver data
588: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
589:   PetscCallA(DMShellCreate(PETSC_COMM_WORLD, dalam, ierr))
590:   PetscCallA(DMShellSetGlobalVector(dalam, x2, ierr))
591:   PetscCallA(DMShellSetMatrix(dalam, Dmat, ierr))

593:   PetscCallA(VecCreate(PETSC_COMM_SELF, x2loc, ierr))
594:   PetscCallA(VecSetSizes(x2loc, nloclam, nloclam, ierr))
595:   PetscCallA(VecSetFromOptions(x2loc, ierr))
596:   PetscCallA(DMShellSetLocalVector(dalam, x2loc, ierr))

598:   PetscCallA(DMSetOptionsPrefix(dalam, 'lambda_', ierr))
599:   PetscCallA(DMSetFromOptions(dalam, ierr))
600: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
601: !  Create field split DA
602: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
603:   PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD, solver%da, ierr))
604:   PetscCallA(DMCompositeAddDM(solver%da, daphi, ierr))
605:   PetscCallA(DMCompositeAddDM(solver%da, dalam, ierr))
606:   PetscCallA(DMSetFromOptions(solver%da, ierr))
607:   PetscCallA(DMSetUp(solver%da, ierr))
608:   PetscCallA(DMCompositeGetGlobalISs(solver%da, isglobal, ierr))
609:   solver%isPhi = isglobal(1)
610:   PetscCallA(PetscObjectReference(solver%isPhi, ierr))
611:   solver%isLambda = isglobal(2)
612:   PetscCallA(PetscObjectReference(solver%isLambda, ierr))

614: !     cache matrices
615:   solver%Amat = Amat
616:   solver%Bmat = Bmat
617:   solver%Cmat = Cmat
618:   solver%Dmat = Dmat

620:   matArray(1) = Amat
621:   matArray(2) = Bmat
622:   matArray(3) = Cmat
623:   matArray(4) = Dmat

625:   PetscCallA(MatCreateNest(PETSC_COMM_WORLD, itwo, isglobal, itwo, isglobal, matArray, KKTmat, ierr))
626:   PetscCallA(DMCompositeRestoreGlobalISs(solver%da, isglobal, ierr))
627:   PetscCallA(MatSetFromOptions(KKTmat, ierr))

629: !  Extract global and local vectors from DMDA; then duplicate for remaining
630: !     vectors that are the same types
631:   PetscCallA(MatCreateVecs(KKTmat, x, PETSC_NULL_VEC, ierr))
632:   PetscCallA(VecDuplicate(x, r, ierr))

634:   PetscCallA(SNESCreate(PETSC_COMM_WORLD, mysnes, ierr))

636:   PetscCallA(SNESSetDM(mysnes, solver%da, ierr))

638:   PetscCallA(SNESSetApplicationContext(mysnes, solver, ierr))

640:   PetscCallA(SNESSetDM(mysnes, solver%da, ierr))

642: !  Set function evaluation routine and vector
643:   PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver, ierr))
644:   if (useobjective .eqv. PETSC_TRUE) then
645:     PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
646:   end if
647:   PetscCallA(SNESSetJacobian(mysnes, KKTmat, KKTmat, FormJacobian, solver, ierr))

649: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
650: !  Customize nonlinear solver; set runtime options
651: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
652: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
653:   PetscCallA(SNESSetFromOptions(mysnes, ierr))

655: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
656: !  Evaluate initial guess; then solve nonlinear system.
657: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
658: !  Note: The user should initialize the vector, x, with the initial guess
659: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
660: !  to employ an initial guess of zero, the user should explicitly set
661: !  this vector to zero by calling VecSet().

663:   PetscCallA(FormInitialGuess(mysnes, x, ierr))
664:   PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr))
665:   PetscCallA(SNESGetIterationNumber(mysnes, its, ierr))
666:   if (solver%rank == 0) then
667:     write (6, 100) its
668:   end if
669: 100 format('Number of SNES iterations = ', i5)

671: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
672: !  Free work space.  All PETSc objects should be destroyed when they
673: !  are no longer needed.
674: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
675:   PetscCallA(MatDestroy(KKTmat, ierr))
676:   PetscCallA(MatDestroy(Amat, ierr))
677:   PetscCallA(MatDestroy(Dmat, ierr))
678:   PetscCallA(MatDestroy(Bmat, ierr))
679:   PetscCallA(MatDestroy(Cmat, ierr))
680:   PetscCallA(MatDestroy(solver%AmatLin, ierr))
681:   PetscCallA(ISDestroy(solver%isPhi, ierr))
682:   PetscCallA(ISDestroy(solver%isLambda, ierr))
683:   PetscCallA(VecDestroy(x, ierr))
684:   PetscCallA(VecDestroy(x2, ierr))
685:   PetscCallA(VecDestroy(x1, ierr))
686:   PetscCallA(VecDestroy(x1loc, ierr))
687:   PetscCallA(VecDestroy(x2loc, ierr))
688:   PetscCallA(VecDestroy(r, ierr))
689:   PetscCallA(SNESDestroy(mysnes, ierr))
690:   PetscCallA(DMDestroy(solver%da, ierr))
691:   PetscCallA(DMDestroy(daphi, ierr))
692:   PetscCallA(DMDestroy(dalam, ierr))

694:   PetscCallA(PetscFinalize(ierr))
695: end

697: !/*TEST
698: !
699: !   build:
700: !      requires: !single !complex
701: !
702: !   test:
703: !      nsize: 4
704: !      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_esteig_ksp_type cg -fieldsplit_phi_pc_gamg_esteig_ksp_max_it 10 -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0.
705: !
706: !   test:
707: !      args: -snes_linesearch_type {{secant cp}separate output} -objective {{false true}shared output}
708: !
709: !   test:
710: !      args: -snes_linesearch_type bt -objective {{false true}separate output}
711: !
712: !TEST*/