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