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*/

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

62:       end module petsc_kkt_solver

64:       module petsc_kkt_solver_interfaces
65:         use petsc_kkt_solver

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

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

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

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

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

139: !  Initialize problem parameters
140:       lambda_max  = 6.81_PETSC_REAL_KIND
141:       lambda_min  = 0.0
142:       solver%lambda = 6.0
143:       ione = 1
144:       nfour = 4
145:       itwo = 2
146:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr)
147:       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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

367:       call PetscFinalize(ierr)
368:       end

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

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

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

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

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

419:       return
420:       end subroutine FormInitialGuess

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

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

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

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

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

474:       return
475:       end subroutine InitialGuessLocal

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

503: !  Declarations for use with local arrays:
504:       Vec::      Xsub(1)
505:      Mat::     Amat
506:       PetscInt       ione

508:       ione = 1

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

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

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

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

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

531:       return
532:       end subroutine FormJacobian

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

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

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

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

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

617:       return
618:       end subroutine FormJacobianLocal

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

644: !  Declarations for use with local arrays:
645:      Vec::              Xsub(2),Fsub(2)
646:       PetscInt               itwo

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

653:       itwo = 2
654:       call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
655:       call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)

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

660: !     do rest of operator (linear)
661:       call MatMult(    solver%Cmat, Xsub(1),      Fsub(2), ierr);CHKERRQ(ierr)
662:       call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
663:       call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr)

665:       call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
666:       call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
667:       return
668:       end subroutine formfunction

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

700:       sc     = solver%lambda
701:       ione   = 1

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

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

724:       call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr)
725:       call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr)

727:       0
728:       return
729:       end subroutine FormFunctionNLTerm

731: !/*TEST
732: !
733: !   build:
734: !      requires: !single !complex
735: !
736: !   test:
737: !      nsize: 4
738: !      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
739: !
740: !TEST*/
```