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: module ex73f90tmodule
 42: #include <petsc/finclude/petscdmda.h>
 43: #include <petsc/finclude/petscdmcomposite.h>
 44: #include <petsc/finclude/petscmat.h>
 45:   use petscdmda
 46:   use petscdmcomposite
 47:   use petscmat
 48:   type ex73f90tmodule_type
 49:     DM::da
 50: !     temp A block stuff
 51:     PetscInt mx, my
 52:     PetscMPIInt rank
 53:     PetscReal lambda
 54: !     Mats
 55:     Mat::Amat, AmatLin, Bmat, CMat, Dmat
 56:     IS::isPhi, isLambda
 57:   end type ex73f90tmodule_type

 59: end module ex73f90tmodule

 61: module ex73f90tmodule_interfaces
 62:   use ex73f90tmodule

 64:   Interface SNESSetApplicationContext
 65:     Subroutine SNESSetApplicationContext(snesIn, ctx, ierr)
 66: #include <petsc/finclude/petscsnes.h>
 67:       use petscsnes
 68:       use ex73f90tmodule
 69:       SNES::    snesIn
 70:       type(ex73f90tmodule_type) ctx
 71:       PetscErrorCode ierr
 72:     End Subroutine
 73:   End Interface SNESSetApplicationContext

 75:   Interface SNESGetApplicationContext
 76:     Subroutine SNESGetApplicationContext(snesIn, ctx, ierr)
 77: #include <petsc/finclude/petscsnes.h>
 78:       use petscsnes
 79:       use ex73f90tmodule
 80:       SNES::     snesIn
 81:       type(ex73f90tmodule_type), pointer :: ctx
 82:       PetscErrorCode ierr
 83:     End Subroutine
 84:   End Interface SNESGetApplicationContext
 85: end module ex73f90tmodule_interfaces

 87: subroutine MyObjective(snes, x, result, ctx, ierr)
 88: #include <petsc/finclude/petsc.h>
 89:   use petsc
 90:   implicit none
 91:   PetscInt ctx
 92:   Vec x, f
 93:   SNES snes
 94:   PetscErrorCode ierr
 95:   PetscScalar result
 96:   PetscReal fnorm

 98:   PetscCall(VecDuplicate(x, f, ierr))
 99:   PetscCall(SNESComputeFunction(snes, x, f, ierr))
100:   PetscCall(VecNorm(f, NORM_2, fnorm, ierr))
101:   result = .5*fnorm*fnorm
102:   PetscCall(VecDestroy(f, ierr))
103: end subroutine MyObjective

105: program main
106: #include <petsc/finclude/petscdmda.h>
107: #include <petsc/finclude/petscsnes.h>
108:   use petscsnes
109:   use ex73f90tmodule
110:   use ex73f90tmodule_interfaces
111:   implicit none
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: !                   Variable declarations
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: !
116: !  Variables:
117: !     mysnes      - nonlinear solver
118: !     x, r        - solution, residual vectors
119: !     J           - Jacobian matrix
120: !     its         - iterations for convergence
121: !     Nx, Ny      - number of preocessors in x- and y- directions
122: !
123:   SNES::       mysnes
124:   Vec::        x, r, x2, x1, x1loc, x2loc
125:   Mat::       Amat, Bmat, Cmat, Dmat, KKTMat, matArray(4)
126: !      Mat::       tmat
127:   DM::       daphi, dalam
128:   IS, pointer ::        isglobal(:)
129:   PetscErrorCode ierr
130:   PetscInt its, N1, N2, i, j, irow, row(1)
131:   PetscInt col(1), low, high, lamlow, lamhigh
132:   PetscBool flg
133:   PetscInt ione, nfour, itwo, nloc, nloclam
134:   PetscReal lambda_max, lambda_min
135:   type(ex73f90tmodule_type) solver
136:   PetscScalar bval(1), cval(1), one
137:   PetscBool useobjective

139: !  Note: Any user-defined Fortran routines (such as FormJacobian)
140: !  MUST be declared as external.
141:   external FormInitialGuess, FormJacobian, FormFunction, MyObjective

143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: !  Initialize program
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:   PetscCallA(PetscInitialize(ierr))
147:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, solver%rank, ierr))

149: !  Initialize problem parameters
150:   lambda_max = 6.81_PETSC_REAL_KIND
151:   lambda_min = 0.0
152:   solver%lambda = 6.0
153:   ione = 1
154:   nfour = 4
155:   itwo = 2
156:   useobjective = PETSC_FALSE
157:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', solver%lambda, flg, ierr))
158:   PetscCheckA(solver%lambda <= lambda_max .and. solver%lambda >= lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')
159:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-objective', useobjective, flg, ierr))

161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: !  Create vector data structures; set function evaluation routine
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

165: !     just get size
166:   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))
167:   PetscCallA(DMSetFromOptions(daphi, ierr))
168:   PetscCallA(DMSetUp(daphi, ierr))
169:   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))
170:   N1 = solver%my*solver%mx
171:   N2 = solver%my
172:   flg = .false.
173:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_constraints', flg, flg, ierr))
174:   if (flg) then
175:     N2 = 0
176:   end if

178:   PetscCallA(DMDestroy(daphi, ierr))

180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: !  Create matrix data structure; set Jacobian evaluation routine
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:   PetscCallA(DMShellCreate(PETSC_COMM_WORLD, daphi, ierr))
184:   PetscCallA(DMSetOptionsPrefix(daphi, 'phi_', ierr))
185:   PetscCallA(DMSetFromOptions(daphi, ierr))

187:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x1, ierr))
188:   PetscCallA(VecSetSizes(x1, PETSC_DECIDE, N1, ierr))
189:   PetscCallA(VecSetFromOptions(x1, ierr))

191:   PetscCallA(VecGetOwnershipRange(x1, low, high, ierr))
192:   nloc = high - low

194:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr))
195:   PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
196:   PetscCallA(MatSetUp(Amat, ierr))

198:   PetscCallA(MatCreate(PETSC_COMM_WORLD, solver%AmatLin, ierr))
199:   PetscCallA(MatSetSizes(solver%AmatLin, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
200:   PetscCallA(MatSetUp(solver%AmatLin, ierr))

202:   PetscCallA(FormJacobianLocal(x1, solver%AmatLin, solver, .false., ierr))
203:   PetscCallA(MatAssemblyBegin(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
204:   PetscCallA(MatAssemblyEnd(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))

206:   PetscCallA(DMShellSetGlobalVector(daphi, x1, ierr))
207:   PetscCallA(DMShellSetMatrix(daphi, Amat, ierr))

209:   PetscCallA(VecCreate(PETSC_COMM_SELF, x1loc, ierr))
210:   PetscCallA(VecSetSizes(x1loc, nloc, nloc, ierr))
211:   PetscCallA(VecSetFromOptions(x1loc, ierr))
212:   PetscCallA(DMShellSetLocalVector(daphi, x1loc, ierr))

214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: !  Create B, C, & D matrices
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Cmat, ierr))
218:   PetscCallA(MatSetSizes(Cmat, PETSC_DECIDE, PETSC_DECIDE, N2, N1, ierr))
219:   PetscCallA(MatSetUp(Cmat, ierr))
220: !      create data for C and B
221:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Bmat, ierr))
222:   PetscCallA(MatSetSizes(Bmat, PETSC_DECIDE, PETSC_DECIDE, N1, N2, ierr))
223:   PetscCallA(MatSetUp(Bmat, ierr))
224: !     create data for D
225:   PetscCallA(MatCreate(PETSC_COMM_WORLD, Dmat, ierr))
226:   PetscCallA(MatSetSizes(Dmat, PETSC_DECIDE, PETSC_DECIDE, N2, N2, ierr))
227:   PetscCallA(MatSetUp(Dmat, ierr))

229:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x2, ierr))
230:   PetscCallA(VecSetSizes(x2, PETSC_DECIDE, N2, ierr))
231:   PetscCallA(VecSetFromOptions(x2, ierr))

233:   PetscCallA(VecGetOwnershipRange(x2, lamlow, lamhigh, ierr))
234:   nloclam = lamhigh - lamlow

236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: !  Set fake B and C
238: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239:   one = 1.0
240:   if (N2 > 0) then
241:     bval(1) = -one/(solver%mx - 2)
242: !     cval = -one/(solver%my*solver%mx)
243:     cval(1) = -one
244:     do 20 irow = low, high - 1
245:       j = irow/solver%mx   ! row in domain
246:       i = mod(irow, solver%mx)
247:       row(1) = irow
248:       col(1) = j
249:       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
250:         !     no op
251:       else
252:         PetscCallA(MatSetValues(Bmat, ione, row, ione, col, bval, INSERT_VALUES, ierr))
253:       end if
254:       row(1) = j
255:       PetscCallA(MatSetValues(Cmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
256: 20    continue
257:       end if
258:       PetscCallA(MatAssemblyBegin(Bmat, MAT_FINAL_ASSEMBLY, ierr))
259:       PetscCallA(MatAssemblyEnd(Bmat, MAT_FINAL_ASSEMBLY, ierr))
260:       PetscCallA(MatAssemblyBegin(Cmat, MAT_FINAL_ASSEMBLY, ierr))
261:       PetscCallA(MatAssemblyEnd(Cmat, MAT_FINAL_ASSEMBLY, ierr))

263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: !  Set D (identity)
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266:       do 30 j = lamlow, lamhigh - 1
267:         row(1) = j
268:         cval(1) = one
269:         PetscCallA(MatSetValues(Dmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
270: 30      continue
271:         PetscCallA(MatAssemblyBegin(Dmat, MAT_FINAL_ASSEMBLY, ierr))
272:         PetscCallA(MatAssemblyEnd(Dmat, MAT_FINAL_ASSEMBLY, ierr))

274: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: !  DM for lambda (dalam) : temp driver for A block, setup A block solver data
276: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277:         PetscCallA(DMShellCreate(PETSC_COMM_WORLD, dalam, ierr))
278:         PetscCallA(DMShellSetGlobalVector(dalam, x2, ierr))
279:         PetscCallA(DMShellSetMatrix(dalam, Dmat, ierr))

281:         PetscCallA(VecCreate(PETSC_COMM_SELF, x2loc, ierr))
282:         PetscCallA(VecSetSizes(x2loc, nloclam, nloclam, ierr))
283:         PetscCallA(VecSetFromOptions(x2loc, ierr))
284:         PetscCallA(DMShellSetLocalVector(dalam, x2loc, ierr))

286:         PetscCallA(DMSetOptionsPrefix(dalam, 'lambda_', ierr))
287:         PetscCallA(DMSetFromOptions(dalam, ierr))
288: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: !  Create field split DA
290: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291:         PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD, solver%da, ierr))
292:         PetscCallA(DMCompositeAddDM(solver%da, daphi, ierr))
293:         PetscCallA(DMCompositeAddDM(solver%da, dalam, ierr))
294:         PetscCallA(DMSetFromOptions(solver%da, ierr))
295:         PetscCallA(DMSetUp(solver%da, ierr))
296:         PetscCallA(DMCompositeGetGlobalISs(solver%da, isglobal, ierr))
297:         solver%isPhi = isglobal(1)
298:         PetscCallA(PetscObjectReference(solver%isPhi, ierr))
299:         solver%isLambda = isglobal(2)
300:         PetscCallA(PetscObjectReference(solver%isLambda, ierr))

302: !     cache matrices
303:         solver%Amat = Amat
304:         solver%Bmat = Bmat
305:         solver%Cmat = Cmat
306:         solver%Dmat = Dmat

308:         matArray(1) = Amat
309:         matArray(2) = Bmat
310:         matArray(3) = Cmat
311:         matArray(4) = Dmat

313:         PetscCallA(MatCreateNest(PETSC_COMM_WORLD, itwo, isglobal, itwo, isglobal, matArray, KKTmat, ierr))
314:         PetscCallA(DMCompositeRestoreGlobalISs(solver%da, isglobal, ierr))
315:         PetscCallA(MatSetFromOptions(KKTmat, ierr))

317: !  Extract global and local vectors from DMDA; then duplicate for remaining
318: !     vectors that are the same types
319:         PetscCallA(MatCreateVecs(KKTmat, x, PETSC_NULL_VEC, ierr))
320:         PetscCallA(VecDuplicate(x, r, ierr))

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

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

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

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

330: !  Set function evaluation routine and vector
331:         PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver, ierr))
332:         if (useobjective .eqv. PETSC_TRUE) then
333:           PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
334:         end if
335:         PetscCallA(SNESSetJacobian(mysnes, KKTmat, KKTmat, FormJacobian, solver, ierr))

337: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
338: !  Customize nonlinear solver; set runtime options
339: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
341:         PetscCallA(SNESSetFromOptions(mysnes, ierr))

343: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344: !  Evaluate initial guess; then solve nonlinear system.
345: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
346: !  Note: The user should initialize the vector, x, with the initial guess
347: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
348: !  to employ an initial guess of zero, the user should explicitly set
349: !  this vector to zero by calling VecSet().

351:         PetscCallA(FormInitialGuess(mysnes, x, ierr))
352:         PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr))
353:         PetscCallA(SNESGetIterationNumber(mysnes, its, ierr))
354:         if (solver%rank == 0) then
355:           write (6, 100) its
356:         end if
357: 100     format('Number of SNES iterations = ', i5)

359: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
360: !  Free work space.  All PETSc objects should be destroyed when they
361: !  are no longer needed.
362: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
363:         PetscCallA(MatDestroy(KKTmat, ierr))
364:         PetscCallA(MatDestroy(Amat, ierr))
365:         PetscCallA(MatDestroy(Dmat, ierr))
366:         PetscCallA(MatDestroy(Bmat, ierr))
367:         PetscCallA(MatDestroy(Cmat, ierr))
368:         PetscCallA(MatDestroy(solver%AmatLin, ierr))
369:         PetscCallA(ISDestroy(solver%isPhi, ierr))
370:         PetscCallA(ISDestroy(solver%isLambda, ierr))
371:         PetscCallA(VecDestroy(x, ierr))
372:         PetscCallA(VecDestroy(x2, ierr))
373:         PetscCallA(VecDestroy(x1, ierr))
374:         PetscCallA(VecDestroy(x1loc, ierr))
375:         PetscCallA(VecDestroy(x2loc, ierr))
376:         PetscCallA(VecDestroy(r, ierr))
377:         PetscCallA(SNESDestroy(mysnes, ierr))
378:         PetscCallA(DMDestroy(solver%da, ierr))
379:         PetscCallA(DMDestroy(daphi, ierr))
380:         PetscCallA(DMDestroy(dalam, ierr))

382:         PetscCallA(PetscFinalize(ierr))
383:       end

385: ! ---------------------------------------------------------------------
386: !
387: !  FormInitialGuess - Forms initial approximation.
388: !
389: !  Input Parameters:
390: !  X - vector
391: !
392: !  Output Parameter:
393: !  X - vector
394: !
395: !  Notes:
396: !  This routine serves as a wrapper for the lower-level routine
397: !  "InitialGuessLocal", where the actual computations are
398: !  done using the standard Fortran style of treating the local
399: !  vector data as a multidimensional array over the local mesh.
400: !  This routine merely handles ghost point scatters and accesses
401: !  the local vector data via VecGetArray() and VecRestoreArray().
402: !
403:       subroutine FormInitialGuess(mysnes, Xnest, ierr)
404: #include <petsc/finclude/petscsnes.h>
405:         use petscsnes
406:         use ex73f90tmodule
407:         use ex73f90tmodule_interfaces
408:         implicit none
409: !  Input/output variables:
410:         SNES::     mysnes
411:         Vec::      Xnest
412:         PetscErrorCode ierr

414: !  Declarations for use with local arrays:
415:         type(ex73f90tmodule_type), pointer:: solver
416:         Vec::      Xsub(2)
417:         PetscInt::  izero, ione, itwo

419:         izero = 0
420:         ione = 1
421:         itwo = 2
422:         ierr = 0
423:         PetscCall(SNESGetApplicationContext(mysnes, solver, ierr))
424:         PetscCall(DMCompositeGetAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

426:         PetscCall(InitialGuessLocal(solver, Xsub(1), ierr))
427:         PetscCall(VecAssemblyBegin(Xsub(1), ierr))
428:         PetscCall(VecAssemblyEnd(Xsub(1), ierr))

430: !     zero out lambda
431:         PetscCall(VecZeroEntries(Xsub(2), ierr))
432:         PetscCall(DMCompositeRestoreAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))

434:       end subroutine FormInitialGuess

436: ! ---------------------------------------------------------------------
437: !
438: !  InitialGuessLocal - Computes initial approximation, called by
439: !  the higher level routine FormInitialGuess().
440: !
441: !  Input Parameter:
442: !  X1 - local vector data
443: !
444: !  Output Parameters:
445: !  x - local vector data
446: !  ierr - error code
447: !
448: !  Notes:
449: !  This routine uses standard Fortran-style computations over a 2-dim array.
450: !
451:       subroutine InitialGuessLocal(solver, X1, ierr)
452: #include <petsc/finclude/petscsys.h>
453:         use petscsys
454:         use ex73f90tmodule
455:         implicit none
456: !  Input/output variables:
457:         type(ex73f90tmodule_type) solver
458:         Vec::      X1
459:         PetscErrorCode ierr

461: !  Local variables:
462:         PetscInt row, i, j, ione, low, high
463:         PetscReal temp1, temp, hx, hy, v
464:         PetscReal one

466: !  Set parameters
467:         ione = 1
468:         ierr = 0
469:         one = 1.0
470:         hx = one/(solver%mx - 1)
471:         hy = one/(solver%my - 1)
472:         temp1 = solver%lambda/(solver%lambda + one) + one

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

476:         do 20 row = low, high - 1
477:           j = row/solver%mx
478:           i = mod(row, solver%mx)
479:           temp = min(j, solver%my - j + 1)*hy
480:           if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
481:             v = 0.0
482:           else
483:             v = temp1*sqrt(min(min(i, solver%mx - i + 1)*hx, temp))
484:           end if
485:           PetscCall(VecSetValues(X1, ione, [row], [v], INSERT_VALUES, ierr))
486: 20        continue

488:           end subroutine InitialGuessLocal

490: ! ---------------------------------------------------------------------
491: !
492: !  FormJacobian - Evaluates Jacobian matrix.
493: !
494: !  Input Parameters:
495: !  dummy     - the SNES context
496: !  x         - input vector
497: !  solver    - solver data
498: !
499: !  Output Parameters:
500: !  jac      - Jacobian matrix
501: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
502: !
503:           subroutine FormJacobian(dummy, X, jac, jac_prec, solver, ierr)
504: #include <petsc/finclude/petscsnes.h>
505:             use petscsnes
506:             use ex73f90tmodule
507:             implicit none
508: !  Input/output variables:
509:             SNES::     dummy
510:             Vec::      X
511:             Mat::     jac, jac_prec
512:             type(ex73f90tmodule_type) solver
513:             PetscErrorCode ierr

515: !  Declarations for use with local arrays:
516:             Vec::      Xsub(1)
517:             Mat::     Amat
518:             PetscInt ione

520:             ione = 1

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

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

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

531:             ! the rest of the matrix is not touched
532:             PetscCall(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
533:             PetscCall(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
534:             if (jac /= jac_prec) then
535:               PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
536:               PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
537:             end if

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

543:           end subroutine FormJacobian

545: ! ---------------------------------------------------------------------
546: !
547: !  FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
548: !  called by the higher level routine FormJacobian().
549: !
550: !  Input Parameters:
551: !  x        - local vector data
552: !
553: !  Output Parameters:
554: !  jac - Jacobian matrix used to compute the preconditioner
555: !  ierr     - error code
556: !
557: !  Notes:
558: !  This routine uses standard Fortran-style computations over a 2-dim array.
559: !
560:           subroutine FormJacobianLocal(X1, jac, solver, add_nl_term, ierr)
561: #include <petsc/finclude/petscmat.h>
562:             use ex73f90tmodule
563:             implicit none
564: !  Input/output variables:
565:             type(ex73f90tmodule_type) solver
566:             Vec::      X1
567:             Mat::     jac
568:             logical add_nl_term
569:             PetscErrorCode ierr

571: !  Local variables:
572:             PetscInt irow, row(1), col(5), i, j
573:             PetscInt ione, ifive, low, high, ii
574:             PetscScalar two, one, hx, hy, hy2inv
575:             PetscScalar hx2inv, sc, v(5)
576:             PetscScalar, pointer :: lx_v(:)

578: !  Set parameters
579:             ione = 1
580:             ifive = 5
581:             one = 1.0
582:             two = 2.0
583:             hx = one/(solver%mx - 1)
584:             hy = one/(solver%my - 1)
585:             sc = solver%lambda
586:             hx2inv = one/(hx*hx)
587:             hy2inv = one/(hy*hy)

589:             PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
590:             PetscCall(VecGetArrayRead(X1, lx_v, ierr))

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

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

627:               end subroutine FormJacobianLocal

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

653: !  Declarations for use with local arrays:
654:                 Vec::              Xsub(2), Fsub(2)
655:                 PetscInt itwo

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

662:                 itwo = 2
663:                 PetscCall(DMCompositeGetAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
664:                 PetscCall(DMCompositeGetAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))

666:                 PetscCall(FormFunctionNLTerm(Xsub(1), Fsub(1), solver, ierr))
667:                 PetscCall(MatMultAdd(solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))

669: !     do rest of operator (linear)
670:                 PetscCall(MatMult(solver%Cmat, Xsub(1), Fsub(2), ierr))
671:                 PetscCall(MatMultAdd(solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
672:                 PetscCall(MatMultAdd(solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))

674:                 PetscCall(DMCompositeRestoreAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
675:                 PetscCall(DMCompositeRestoreAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
676:               end subroutine formfunction

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

707:                 sc = solver%lambda
708:                 ione = 1

710:                 PetscCall(VecGetArrayRead(X1, lx_v, ierr))
711:                 PetscCall(VecGetOwnershipRange(X1, low, high, ierr))

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

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

731:                   PetscCall(VecAssemblyBegin(F1, ierr))
732:                   PetscCall(VecAssemblyEnd(F1, ierr))

734:                   ierr = 0
735:                   end subroutine FormFunctionNLTerm

737: !/*TEST
738: !
739: !   build:
740: !      requires: !single !complex
741: !
742: !   test:
743: !      nsize: 4
744: !      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.
745: !
746: !   test:
747: !      args: -snes_linesearch_type {{secant cp}separate output} -objective {{false true}shared output}
748: !
749: !   test:
750: !      args: -snes_linesearch_type bt -objective {{false true}separate output}
751: !
752: !TEST*/