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: contains
61: subroutine MyObjective(snes, x, result, ctx, ierr)
62: PetscInt ctx
63: Vec x, f
64: SNES snes
65: PetscErrorCode ierr
66: PetscScalar result
67: PetscReal fnorm
69: PetscCall(VecDuplicate(x, f, ierr))
70: PetscCall(SNESComputeFunction(snes, x, f, ierr))
71: PetscCall(VecNorm(f, NORM_2, fnorm, ierr))
72: result = .5*fnorm*fnorm
73: PetscCall(VecDestroy(f, ierr))
74: end subroutine MyObjective
76: ! ---------------------------------------------------------------------
77: !
78: ! FormInitialGuess - Forms initial approximation.
79: !
80: ! Input Parameters:
81: ! X - vector
82: !
83: ! Output Parameter:
84: ! X - vector
85: !
86: ! Notes:
87: ! This routine serves as a wrapper for the lower-level routine
88: ! "InitialGuessLocal", where the actual computations are
89: ! done using the standard Fortran style of treating the local
90: ! vector data as a multidimensional array over the local mesh.
91: ! This routine merely handles ghost point scatters and accesses
92: ! the local vector data via VecGetArray() and VecRestoreArray().
93: !
94: subroutine FormInitialGuess(mysnes, Xnest, ierr)
95: ! Input/output variables:
96: SNES:: mysnes
97: Vec:: Xnest
98: PetscErrorCode ierr
100: ! Declarations for use with local arrays:
101: type(ex73f90tmodule_type), pointer:: solver
102: Vec:: Xsub(2)
103: PetscInt:: izero, ione, itwo
105: izero = 0
106: ione = 1
107: itwo = 2
108: ierr = 0
109: PetscCall(SNESGetApplicationContext(mysnes, solver, ierr))
110: PetscCall(DMCompositeGetAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
112: PetscCall(InitialGuessLocal(solver, Xsub(1), ierr))
113: PetscCall(VecAssemblyBegin(Xsub(1), ierr))
114: PetscCall(VecAssemblyEnd(Xsub(1), ierr))
116: ! zero out lambda
117: PetscCall(VecZeroEntries(Xsub(2), ierr))
118: PetscCall(DMCompositeRestoreAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
120: end subroutine FormInitialGuess
122: ! ---------------------------------------------------------------------
123: !
124: ! InitialGuessLocal - Computes initial approximation, called by
125: ! the higher level routine FormInitialGuess().
126: !
127: ! Input Parameter:
128: ! X1 - local vector data
129: !
130: ! Output Parameters:
131: ! x - local vector data
132: ! ierr - error code
133: !
134: ! Notes:
135: ! This routine uses standard Fortran-style computations over a 2-dim array.
136: !
137: subroutine InitialGuessLocal(solver, X1, ierr)
138: ! Input/output variables:
139: type(ex73f90tmodule_type) solver
140: Vec:: X1
141: PetscErrorCode ierr
143: ! Local variables:
144: PetscInt row, i, j, ione, low, high
145: PetscReal temp1, temp, hx, hy, v
146: PetscReal one
148: ! Set parameters
149: ione = 1
150: ierr = 0
151: one = 1.0
152: hx = one/(solver%mx - 1)
153: hy = one/(solver%my - 1)
154: temp1 = solver%lambda/(solver%lambda + one) + one
156: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
158: do row = low, high - 1
159: j = row/solver%mx
160: i = mod(row, solver%mx)
161: temp = min(j, solver%my - j + 1)*hy
162: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
163: v = 0.0
164: else
165: v = temp1*sqrt(min(min(i, solver%mx - i + 1)*hx, temp))
166: end if
167: PetscCall(VecSetValues(X1, ione, [row], [v], INSERT_VALUES, ierr))
168: end do
170: end subroutine InitialGuessLocal
172: ! ---------------------------------------------------------------------
173: !
174: ! FormJacobian - Evaluates Jacobian matrix.
175: !
176: ! Input Parameters:
177: ! dummy - the SNES context
178: ! x - input vector
179: ! solver - solver data
180: !
181: ! Output Parameters:
182: ! jac - Jacobian matrix
183: ! jac_prec - optionally different matrix used to construct the preconditioner (not used here)
184: !
185: subroutine FormJacobian(dummy, X, jac, jac_prec, solver, ierr)
186: ! Input/output variables:
187: SNES:: dummy
188: Vec:: X
189: Mat:: jac, jac_prec
190: type(ex73f90tmodule_type) solver
191: PetscErrorCode ierr
193: ! Declarations for use with local arrays:
194: Vec:: Xsub(1)
195: Mat:: Amat
196: PetscInt ione
198: ione = 1
200: PetscCall(DMCompositeGetAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
202: ! Compute entries for the locally owned part of the Jacobian preconditioner.
203: PetscCall(MatCreateSubMatrix(jac_prec, solver%isPhi, solver%isPhi, MAT_INITIAL_MATRIX, Amat, ierr))
205: PetscCall(FormJacobianLocal(Xsub(1), Amat, solver, .true., ierr))
206: PetscCall(MatDestroy(Amat, ierr)) ! discard our reference
207: PetscCall(DMCompositeRestoreAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
209: ! the rest of the matrix is not touched
210: PetscCall(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
211: PetscCall(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
212: if (jac /= jac_prec) then
213: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
214: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
215: end if
217: ! Tell the matrix we will never add a new nonzero location to the
218: ! matrix. If we do it will generate an error.
219: PetscCall(MatSetOption(jac_prec, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
221: end subroutine FormJacobian
223: ! ---------------------------------------------------------------------
224: !
225: ! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
226: ! called by the higher level routine FormJacobian().
227: !
228: ! Input Parameters:
229: ! x - local vector data
230: !
231: ! Output Parameters:
232: ! jac - Jacobian matrix used to compute the preconditioner
233: ! ierr - error code
234: !
235: ! Notes:
236: ! This routine uses standard Fortran-style computations over a 2-dim array.
237: !
238: subroutine FormJacobianLocal(X1, jac, solver, add_nl_term, ierr)
239: ! Input/output variables:
240: type(ex73f90tmodule_type) solver
241: Vec:: X1
242: Mat:: jac
243: logical add_nl_term
244: PetscErrorCode ierr
246: ! Local variables:
247: PetscInt irow, row(1), col(5), i, j
248: PetscInt ione, ifive, low, high, ii
249: PetscScalar two, one, hx, hy, hy2inv
250: PetscScalar hx2inv, sc, v(5)
251: PetscScalar, pointer :: lx_v(:)
253: ! Set parameters
254: ione = 1
255: ifive = 5
256: one = 1.0
257: two = 2.0
258: hx = one/(solver%mx - 1)
259: hy = one/(solver%my - 1)
260: sc = solver%lambda
261: hx2inv = one/(hx*hx)
262: hy2inv = one/(hy*hy)
264: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
265: PetscCall(VecGetArrayRead(X1, lx_v, ierr))
267: ii = 0
268: do irow = low, high - 1
269: j = irow/solver%mx
270: i = mod(irow, solver%mx)
271: ii = ii + 1 ! one based local index
272: ! boundary points
273: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
274: col(1) = irow
275: row(1) = irow
276: v(1) = one
277: PetscCall(MatSetValues(jac, ione, row, ione, col, v, INSERT_VALUES, ierr))
278: ! interior grid points
279: else
280: v(1) = -hy2inv
281: if (j - 1 == 0) v(1) = 0.0
282: v(2) = -hx2inv
283: if (i - 1 == 0) v(2) = 0.0
284: v(3) = two*(hx2inv + hy2inv)
285: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
286: v(4) = -hx2inv
287: if (i + 1 == solver%mx - 1) v(4) = 0.0
288: v(5) = -hy2inv
289: if (j + 1 == solver%my - 1) v(5) = 0.0
290: col(1) = irow - solver%mx
291: col(2) = irow - 1
292: col(3) = irow
293: col(4) = irow + 1
294: col(5) = irow + solver%mx
295: row(1) = irow
296: PetscCall(MatSetValues(jac, ione, row, ifive, col, v, INSERT_VALUES, ierr))
297: end if
298: end do
300: PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))
302: end subroutine FormJacobianLocal
304: ! ---------------------------------------------------------------------
305: !
306: ! FormFunction - Evaluates nonlinear function, F(x).
307: !
308: ! Input Parameters:
309: ! snes - the SNES context
310: ! X - input vector
311: ! dummy - optional user-defined context, as set by SNESSetFunction()
312: ! (not used here)
313: !
314: ! Output Parameter:
315: ! F - function vector
316: !
317: subroutine FormFunction(snesIn, X, F, solver, ierr)
318: ! Input/output variables:
319: SNES:: snesIn
320: Vec:: X, F
321: PetscErrorCode ierr
322: type(ex73f90tmodule_type) solver
324: ! Declarations for use with local arrays:
325: Vec:: Xsub(2), Fsub(2)
326: PetscInt itwo
328: ! Scatter ghost points to local vector, using the 2-step process
329: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
330: ! By placing code between these two statements, computations can
331: ! be done while messages are in transition.
333: itwo = 2
334: PetscCall(DMCompositeGetAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
335: PetscCall(DMCompositeGetAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
337: PetscCall(FormFunctionNLTerm(Xsub(1), Fsub(1), solver, ierr))
338: PetscCall(MatMultAdd(solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
340: ! do rest of operator (linear)
341: PetscCall(MatMult(solver%Cmat, Xsub(1), Fsub(2), ierr))
342: PetscCall(MatMultAdd(solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
343: PetscCall(MatMultAdd(solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
345: PetscCall(DMCompositeRestoreAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
346: PetscCall(DMCompositeRestoreAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
347: end subroutine formfunction
349: ! ---------------------------------------------------------------------
350: !
351: ! FormFunctionNLTerm - Computes nonlinear function, called by
352: ! the higher level routine FormFunction().
353: !
354: ! Input Parameter:
355: ! x - local vector data
356: !
357: ! Output Parameters:
358: ! f - local vector data, f(x)
359: ! ierr - error code
360: !
361: ! Notes:
362: ! This routine uses standard Fortran-style computations over a 2-dim array.
363: !
364: subroutine FormFunctionNLTerm(X1, F1, solver, ierr)
365: ! Input/output variables:
366: type(ex73f90tmodule_type) solver
367: Vec:: X1, F1
368: PetscErrorCode ierr
369: ! Local variables:
370: PetscScalar sc
371: PetscScalar u, v(1)
372: PetscInt i, j, low, high, ii, ione, irow, row(1)
373: PetscScalar, pointer :: lx_v(:)
375: sc = solver%lambda
376: ione = 1
378: PetscCall(VecGetArrayRead(X1, lx_v, ierr))
379: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
381: ! Compute function over the locally owned part of the grid
382: ii = 0
383: do irow = low, high - 1
384: j = irow/solver%mx
385: i = mod(irow, solver%mx)
386: ii = ii + 1 ! one based local index
387: row(1) = irow
388: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
389: v(1) = 0.0
390: else
391: u = lx_v(ii)
392: v(1) = -sc*exp(u)
393: end if
394: PetscCall(VecSetValues(F1, ione, row, v, INSERT_VALUES, ierr))
395: end do
397: PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))
399: PetscCall(VecAssemblyBegin(F1, ierr))
400: PetscCall(VecAssemblyEnd(F1, ierr))
402: ierr = 0
403: end subroutine FormFunctionNLTerm
405: end module ex73f90tmodule
407: program main
408: use petscsnes
409: use ex73f90tmodule
410: implicit none
411: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
412: ! Variable declarations
413: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414: !
415: ! Variables:
416: ! mysnes - nonlinear solver
417: ! x, r - solution, residual vectors
418: ! J - Jacobian matrix
419: ! its - iterations for convergence
420: ! Nx, Ny - number of preocessors in x- and y- directions
421: !
422: SNES:: mysnes
423: Vec:: x, r, x2, x1, x1loc, x2loc
424: Mat:: Amat, Bmat, Cmat, Dmat, KKTMat, matArray(4)
425: ! Mat:: tmat
426: DM:: daphi, dalam
427: IS, pointer :: isglobal(:)
428: PetscErrorCode ierr
429: PetscInt its, N1, N2, i, j, irow, row(1)
430: PetscInt col(1), low, high, lamlow, lamhigh
431: PetscBool flg
432: PetscInt ione, nfour, itwo, nloc, nloclam
433: PetscReal lambda_max, lambda_min
434: type(ex73f90tmodule_type) solver
435: PetscScalar bval(1), cval(1), one
436: PetscBool useobjective
438: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439: ! Initialize program
440: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
441: PetscCallA(PetscInitialize(ierr))
442: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, solver%rank, ierr))
444: ! Initialize problem parameters
445: lambda_max = 6.81_PETSC_REAL_KIND
446: lambda_min = 0.0
447: solver%lambda = 6.0
448: ione = 1
449: nfour = 4
450: itwo = 2
451: useobjective = PETSC_FALSE
452: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', solver%lambda, flg, ierr))
453: PetscCheckA(solver%lambda <= lambda_max .and. solver%lambda >= lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')
454: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-objective', useobjective, flg, ierr))
456: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457: ! Create vector data structures; set function evaluation routine
458: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: ! just get size
461: 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))
462: PetscCallA(DMSetFromOptions(daphi, ierr))
463: PetscCallA(DMSetUp(daphi, ierr))
464: 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))
465: N1 = solver%my*solver%mx
466: N2 = solver%my
467: flg = .false.
468: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_constraints', flg, flg, ierr))
469: if (flg) then
470: N2 = 0
471: end if
473: PetscCallA(DMDestroy(daphi, ierr))
475: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
476: ! Create matrix data structure; set Jacobian evaluation routine
477: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
478: PetscCallA(DMShellCreate(PETSC_COMM_WORLD, daphi, ierr))
479: PetscCallA(DMSetOptionsPrefix(daphi, 'phi_', ierr))
480: PetscCallA(DMSetFromOptions(daphi, ierr))
482: PetscCallA(VecCreate(PETSC_COMM_WORLD, x1, ierr))
483: PetscCallA(VecSetSizes(x1, PETSC_DECIDE, N1, ierr))
484: PetscCallA(VecSetFromOptions(x1, ierr))
486: PetscCallA(VecGetOwnershipRange(x1, low, high, ierr))
487: nloc = high - low
489: PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr))
490: PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
491: PetscCallA(MatSetUp(Amat, ierr))
493: PetscCallA(MatCreate(PETSC_COMM_WORLD, solver%AmatLin, ierr))
494: PetscCallA(MatSetSizes(solver%AmatLin, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
495: PetscCallA(MatSetUp(solver%AmatLin, ierr))
497: PetscCallA(FormJacobianLocal(x1, solver%AmatLin, solver, .false., ierr))
498: PetscCallA(MatAssemblyBegin(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
499: PetscCallA(MatAssemblyEnd(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
501: PetscCallA(DMShellSetGlobalVector(daphi, x1, ierr))
502: PetscCallA(DMShellSetMatrix(daphi, Amat, ierr))
504: PetscCallA(VecCreate(PETSC_COMM_SELF, x1loc, ierr))
505: PetscCallA(VecSetSizes(x1loc, nloc, nloc, ierr))
506: PetscCallA(VecSetFromOptions(x1loc, ierr))
507: PetscCallA(DMShellSetLocalVector(daphi, x1loc, ierr))
509: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510: ! Create B, C, & D matrices
511: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
512: PetscCallA(MatCreate(PETSC_COMM_WORLD, Cmat, ierr))
513: PetscCallA(MatSetSizes(Cmat, PETSC_DECIDE, PETSC_DECIDE, N2, N1, ierr))
514: PetscCallA(MatSetUp(Cmat, ierr))
515: ! create data for C and B
516: PetscCallA(MatCreate(PETSC_COMM_WORLD, Bmat, ierr))
517: PetscCallA(MatSetSizes(Bmat, PETSC_DECIDE, PETSC_DECIDE, N1, N2, ierr))
518: PetscCallA(MatSetUp(Bmat, ierr))
519: ! create data for D
520: PetscCallA(MatCreate(PETSC_COMM_WORLD, Dmat, ierr))
521: PetscCallA(MatSetSizes(Dmat, PETSC_DECIDE, PETSC_DECIDE, N2, N2, ierr))
522: PetscCallA(MatSetUp(Dmat, ierr))
524: PetscCallA(VecCreate(PETSC_COMM_WORLD, x2, ierr))
525: PetscCallA(VecSetSizes(x2, PETSC_DECIDE, N2, ierr))
526: PetscCallA(VecSetFromOptions(x2, ierr))
528: PetscCallA(VecGetOwnershipRange(x2, lamlow, lamhigh, ierr))
529: nloclam = lamhigh - lamlow
531: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
532: ! Set fake B and C
533: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
534: one = 1.0
535: if (N2 > 0) then
536: bval(1) = -one/(solver%mx - 2)
537: ! cval = -one/(solver%my*solver%mx)
538: cval(1) = -one
539: do irow = low, high - 1
540: j = irow/solver%mx ! row in domain
541: i = mod(irow, solver%mx)
542: row(1) = irow
543: col(1) = j
544: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
545: ! no op
546: else
547: PetscCallA(MatSetValues(Bmat, ione, row, ione, col, bval, INSERT_VALUES, ierr))
548: end if
549: row(1) = j
550: PetscCallA(MatSetValues(Cmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
551: end do
552: end if
553: PetscCallA(MatAssemblyBegin(Bmat, MAT_FINAL_ASSEMBLY, ierr))
554: PetscCallA(MatAssemblyEnd(Bmat, MAT_FINAL_ASSEMBLY, ierr))
555: PetscCallA(MatAssemblyBegin(Cmat, MAT_FINAL_ASSEMBLY, ierr))
556: PetscCallA(MatAssemblyEnd(Cmat, MAT_FINAL_ASSEMBLY, ierr))
558: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
559: ! Set D (identity)
560: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
561: do j = lamlow, lamhigh - 1
562: row(1) = j
563: cval(1) = one
564: PetscCallA(MatSetValues(Dmat, ione, row, ione, row, cval, INSERT_VALUES, ierr))
565: end do
566: PetscCallA(MatAssemblyBegin(Dmat, MAT_FINAL_ASSEMBLY, ierr))
567: PetscCallA(MatAssemblyEnd(Dmat, MAT_FINAL_ASSEMBLY, ierr))
569: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
570: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
571: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
572: PetscCallA(DMShellCreate(PETSC_COMM_WORLD, dalam, ierr))
573: PetscCallA(DMShellSetGlobalVector(dalam, x2, ierr))
574: PetscCallA(DMShellSetMatrix(dalam, Dmat, ierr))
576: PetscCallA(VecCreate(PETSC_COMM_SELF, x2loc, ierr))
577: PetscCallA(VecSetSizes(x2loc, nloclam, nloclam, ierr))
578: PetscCallA(VecSetFromOptions(x2loc, ierr))
579: PetscCallA(DMShellSetLocalVector(dalam, x2loc, ierr))
581: PetscCallA(DMSetOptionsPrefix(dalam, 'lambda_', ierr))
582: PetscCallA(DMSetFromOptions(dalam, ierr))
583: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
584: ! Create field split DA
585: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
586: PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD, solver%da, ierr))
587: PetscCallA(DMCompositeAddDM(solver%da, daphi, ierr))
588: PetscCallA(DMCompositeAddDM(solver%da, dalam, ierr))
589: PetscCallA(DMSetFromOptions(solver%da, ierr))
590: PetscCallA(DMSetUp(solver%da, ierr))
591: PetscCallA(DMCompositeGetGlobalISs(solver%da, isglobal, ierr))
592: solver%isPhi = isglobal(1)
593: PetscCallA(PetscObjectReference(solver%isPhi, ierr))
594: solver%isLambda = isglobal(2)
595: PetscCallA(PetscObjectReference(solver%isLambda, ierr))
597: ! cache matrices
598: solver%Amat = Amat
599: solver%Bmat = Bmat
600: solver%Cmat = Cmat
601: solver%Dmat = Dmat
603: matArray(1) = Amat
604: matArray(2) = Bmat
605: matArray(3) = Cmat
606: matArray(4) = Dmat
608: PetscCallA(MatCreateNest(PETSC_COMM_WORLD, itwo, isglobal, itwo, isglobal, matArray, KKTmat, ierr))
609: PetscCallA(DMCompositeRestoreGlobalISs(solver%da, isglobal, ierr))
610: PetscCallA(MatSetFromOptions(KKTmat, ierr))
612: ! Extract global and local vectors from DMDA; then duplicate for remaining
613: ! vectors that are the same types
614: PetscCallA(MatCreateVecs(KKTmat, x, PETSC_NULL_VEC, ierr))
615: PetscCallA(VecDuplicate(x, r, ierr))
617: PetscCallA(SNESCreate(PETSC_COMM_WORLD, mysnes, ierr))
619: PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
621: PetscCallA(SNESSetApplicationContext(mysnes, solver, ierr))
623: PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
625: ! Set function evaluation routine and vector
626: PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver, ierr))
627: if (useobjective .eqv. PETSC_TRUE) then
628: PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
629: end if
630: PetscCallA(SNESSetJacobian(mysnes, KKTmat, KKTmat, FormJacobian, solver, ierr))
632: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
633: ! Customize nonlinear solver; set runtime options
634: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
635: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
636: PetscCallA(SNESSetFromOptions(mysnes, ierr))
638: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
639: ! Evaluate initial guess; then solve nonlinear system.
640: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
641: ! Note: The user should initialize the vector, x, with the initial guess
642: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
643: ! to employ an initial guess of zero, the user should explicitly set
644: ! this vector to zero by calling VecSet().
646: PetscCallA(FormInitialGuess(mysnes, x, ierr))
647: PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr))
648: PetscCallA(SNESGetIterationNumber(mysnes, its, ierr))
649: if (solver%rank == 0) then
650: write (6, 100) its
651: end if
652: 100 format('Number of SNES iterations = ', i5)
654: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
655: ! Free work space. All PETSc objects should be destroyed when they
656: ! are no longer needed.
657: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
658: PetscCallA(MatDestroy(KKTmat, ierr))
659: PetscCallA(MatDestroy(Amat, ierr))
660: PetscCallA(MatDestroy(Dmat, ierr))
661: PetscCallA(MatDestroy(Bmat, ierr))
662: PetscCallA(MatDestroy(Cmat, ierr))
663: PetscCallA(MatDestroy(solver%AmatLin, ierr))
664: PetscCallA(ISDestroy(solver%isPhi, ierr))
665: PetscCallA(ISDestroy(solver%isLambda, ierr))
666: PetscCallA(VecDestroy(x, ierr))
667: PetscCallA(VecDestroy(x2, ierr))
668: PetscCallA(VecDestroy(x1, ierr))
669: PetscCallA(VecDestroy(x1loc, ierr))
670: PetscCallA(VecDestroy(x2loc, ierr))
671: PetscCallA(VecDestroy(r, ierr))
672: PetscCallA(SNESDestroy(mysnes, ierr))
673: PetscCallA(DMDestroy(solver%da, ierr))
674: PetscCallA(DMDestroy(daphi, ierr))
675: PetscCallA(DMDestroy(dalam, ierr))
677: PetscCallA(PetscFinalize(ierr))
678: end
680: !/*TEST
681: !
682: ! build:
683: ! requires: !single !complex
684: !
685: ! test:
686: ! nsize: 4
687: ! 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.
688: !
689: ! test:
690: ! args: -snes_linesearch_type {{secant cp}separate output} -objective {{false true}shared output}
691: !
692: ! test:
693: ! args: -snes_linesearch_type bt -objective {{false true}separate output}
694: !
695: !TEST*/