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