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, intent(out) :: ierr
99: ! Declarations for use with local arrays:
100: type(ex73f90tmodule_type), pointer:: solver
101: Vec :: Xsub(2)
102: PetscCall(SNESGetApplicationContext(mysnes, solver, ierr))
103: PetscCall(DMCompositeGetAccessArray(solver%da, Xnest, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
105: PetscCall(InitialGuessLocal(solver, Xsub(1), ierr))
106: PetscCall(VecAssemblyBegin(Xsub(1), ierr))
107: PetscCall(VecAssemblyEnd(Xsub(1), ierr))
109: ! zero out lambda
110: PetscCall(VecZeroEntries(Xsub(2), ierr))
111: PetscCall(DMCompositeRestoreAccessArray(solver%da, Xnest, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
113: end subroutine FormInitialGuess
115: ! ---------------------------------------------------------------------
116: !
117: ! InitialGuessLocal - Computes initial approximation, called by
118: ! the higher level routine FormInitialGuess().
119: !
120: ! Input Parameter:
121: ! X1 - local vector data
122: !
123: ! Output Parameters:
124: ! x - local vector data
125: ! ierr - error code
126: !
127: ! Notes:
128: ! This routine uses standard Fortran-style computations over a 2-dim array.
129: !
130: subroutine InitialGuessLocal(solver, X1, ierr)
131: ! Input/output variables:
132: type(ex73f90tmodule_type), intent(in) :: solver
133: Vec :: X1
134: PetscErrorCode, intent(out) :: ierr
135: ! Local variables:
136: PetscInt row, i, j, low, high
137: PetscReal temp1, temp, hx, hy, v
139: ! Set parameters
140: hx = 1._PETSC_REAL_KIND/(solver%mx - 1)
141: hy = 1._PETSC_REAL_KIND/(solver%my - 1)
142: temp1 = solver%lambda/(solver%lambda + 1._PETSC_REAL_KIND) + 1._PETSC_REAL_KIND
144: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
146: do row = low, high - 1
147: j = row/solver%mx
148: i = mod(row, solver%mx)
149: temp = min(j, solver%my - j + 1)*hy
150: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
151: v = 0.0
152: else
153: v = temp1*sqrt(min(min(i, solver%mx - i + 1)*hx, temp))
154: end if
155: PetscCall(VecSetValues(X1, 1_PETSC_INT_KIND, [row], [v], INSERT_VALUES, ierr))
156: end do
158: end subroutine InitialGuessLocal
160: ! ---------------------------------------------------------------------
161: !
162: ! FormJacobian - Evaluates Jacobian matrix.
163: !
164: ! Input Parameters:
165: ! dummy - the SNES context
166: ! x - input vector
167: ! solver - solver data
168: !
169: ! Output Parameters:
170: ! jac - Jacobian matrix
171: ! jac_prec - optionally different matrix used to construct the preconditioner (not used here)
172: !
173: subroutine FormJacobian(dummy, X, jac, jac_prec, solver, ierr)
174: ! Input/output variables:
175: SNES :: dummy
176: Vec :: X
177: Mat :: jac, jac_prec
178: type(ex73f90tmodule_type) solver
179: PetscErrorCode, intent(out) :: ierr
180: ! Declarations for use with local arrays:
181: Vec :: Xsub(1)
182: Mat :: Amat
184: PetscCall(DMCompositeGetAccessArray(solver%da, X, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
186: ! Compute entries for the locally owned part of the Jacobian preconditioner.
187: PetscCall(MatCreateSubMatrix(jac_prec, solver%isPhi, solver%isPhi, MAT_INITIAL_MATRIX, Amat, ierr))
189: PetscCall(FormJacobianLocal(Xsub(1), Amat, solver, .true., ierr))
190: PetscCall(MatDestroy(Amat, ierr)) ! discard our reference
191: PetscCall(DMCompositeRestoreAccessArray(solver%da, X, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
193: ! the rest of the matrix is not touched
194: PetscCall(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
195: PetscCall(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
196: if (jac /= jac_prec) then
197: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
198: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
199: end if
201: ! Tell the matrix we will never add a new nonzero location to the
202: ! matrix. If we do it will generate an error.
203: PetscCall(MatSetOption(jac_prec, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
205: end subroutine FormJacobian
207: ! ---------------------------------------------------------------------
208: !
209: ! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
210: ! called by the higher level routine FormJacobian().
211: !
212: ! Input Parameters:
213: ! x - local vector data
214: !
215: ! Output Parameters:
216: ! jac - Jacobian matrix used to compute the preconditioner
217: ! ierr - error code
218: !
219: ! Notes:
220: ! This routine uses standard Fortran-style computations over a 2-dim array.
221: !
222: subroutine FormJacobianLocal(X1, jac, solver, add_nl_term, ierr)
223: ! Input/output variables:
224: type(ex73f90tmodule_type) solver
225: Vec :: X1
226: Mat :: jac
227: logical add_nl_term
228: PetscErrorCode ierr
229: ! Local variables:
230: PetscInt irow, row(1), col(5), i, j
231: PetscInt low, high, ii
232: PetscScalar, parameter :: one = 1.0, two = 2.0
233: PetscScalar hx, hy, hy2inv, hx2inv, sc, v(5)
234: PetscScalar, pointer :: lx_v(:)
236: ! Set parameters
237: hx = one/(solver%mx - 1)
238: hy = one/(solver%my - 1)
239: sc = solver%lambda
240: hx2inv = one/(hx*hx)
241: hy2inv = one/(hy*hy)
243: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
244: PetscCall(VecGetArrayRead(X1, lx_v, ierr))
246: ii = 0
247: do irow = low, high - 1
248: j = irow/solver%mx
249: i = mod(irow, solver%mx)
250: ii = ii + 1 ! one based local index
251: ! boundary points
252: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
253: col(1) = irow
254: row(1) = irow
255: v(1) = one
256: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
257: ! interior grid points
258: else
259: v(1) = -hy2inv
260: if (j - 1 == 0) v(1) = 0.0
261: v(2) = -hx2inv
262: if (i - 1 == 0) v(2) = 0.0
263: v(3) = two*(hx2inv + hy2inv)
264: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
265: v(4) = -hx2inv
266: if (i + 1 == solver%mx - 1) v(4) = 0.0
267: v(5) = -hy2inv
268: if (j + 1 == solver%my - 1) v(5) = 0.0
269: col(1) = irow - solver%mx
270: col(2) = irow - 1
271: col(3) = irow
272: col(4) = irow + 1
273: col(5) = irow + solver%mx
274: row(1) = irow
275: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, row, 5_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
276: end if
277: end do
279: PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))
281: end subroutine FormJacobianLocal
283: ! ---------------------------------------------------------------------
284: !
285: ! FormFunction - Evaluates nonlinear function, F(x).
286: !
287: ! Input Parameters:
288: ! snes - the SNES context
289: ! X - input vector
290: ! dummy - optional user-defined context, as set by SNESSetFunction()
291: ! (not used here)
292: !
293: ! Output Parameter:
294: ! F - function vector
295: !
296: subroutine FormFunction(snesIn, X, F, solver, ierr)
297: ! Input/output variables:
298: SNES :: snesIn
299: Vec :: X, F
300: PetscErrorCode, intent(out) :: ierr
301: type(ex73f90tmodule_type) solver
303: ! Declarations for use with local arrays:
304: Vec :: Xsub(2), Fsub(2)
306: ! Scatter ghost points to local vector, using the 2-step process
307: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
308: ! By placing code between these two statements, computations can
309: ! be done while messages are in transition.
311: PetscCall(DMCompositeGetAccessArray(solver%da, X, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
312: PetscCall(DMCompositeGetAccessArray(solver%da, F, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
314: PetscCall(FormFunctionNLTerm(Xsub(1), Fsub(1), solver, ierr))
315: PetscCall(MatMultAdd(solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
317: ! do rest of operator (linear)
318: PetscCall(MatMult(solver%Cmat, Xsub(1), Fsub(2), ierr))
319: PetscCall(MatMultAdd(solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
320: PetscCall(MatMultAdd(solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
322: PetscCall(DMCompositeRestoreAccessArray(solver%da, X, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
323: PetscCall(DMCompositeRestoreAccessArray(solver%da, F, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
324: end subroutine formfunction
326: ! ---------------------------------------------------------------------
327: !
328: ! FormFunctionNLTerm - Computes nonlinear function, called by
329: ! the higher level routine FormFunction().
330: !
331: ! Input Parameter:
332: ! x - local vector data
333: !
334: ! Output Parameters:
335: ! f - local vector data, f(x)
336: ! ierr - error code
337: !
338: ! Notes:
339: ! This routine uses standard Fortran-style computations over a 2-dim array.
340: !
341: subroutine FormFunctionNLTerm(X1, F1, solver, ierr)
342: ! Input/output variables:
343: type(ex73f90tmodule_type) solver
344: Vec:: X1, F1
345: PetscErrorCode ierr
346: ! Local variables:
347: PetscScalar sc
348: PetscScalar u, v(1)
349: PetscInt i, j, low, high, ii, irow, row(1)
350: PetscScalar, pointer :: lx_v(:)
352: sc = solver%lambda
354: PetscCall(VecGetArrayRead(X1, lx_v, ierr))
355: PetscCall(VecGetOwnershipRange(X1, low, high, ierr))
357: ! Compute function over the locally owned part of the grid
358: ii = 0
359: do irow = low, high - 1
360: j = irow/solver%mx
361: i = mod(irow, solver%mx)
362: ii = ii + 1 ! one based local index
363: row(1) = irow
364: if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
365: v(1) = 0.0
366: else
367: u = lx_v(ii)
368: v(1) = -sc*exp(u)
369: end if
370: PetscCall(VecSetValues(F1, 1_PETSC_INT_KIND, row, v, INSERT_VALUES, ierr))
371: end do
373: PetscCall(VecRestoreArrayRead(X1, lx_v, ierr))
375: PetscCall(VecAssemblyBegin(F1, ierr))
376: PetscCall(VecAssemblyEnd(F1, ierr))
378: ierr = 0
379: end subroutine FormFunctionNLTerm
381: end module ex73f90tmodule
383: program main
384: use petscsnes
385: use ex73f90tmodule
386: implicit none
387: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: ! Variable declarations
389: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
390: !
391: ! Variables:
392: ! mysnes - nonlinear solver
393: ! x, r - solution, residual vectors
394: ! J - Jacobian matrix
395: ! its - iterations for convergence
396: ! Nx, Ny - number of preocessors in x- and y- directions
397: !
398: SNES :: mysnes
399: Vec :: x, r, x2, x1, x1loc, x2loc
400: Mat :: Amat, Bmat, Cmat, Dmat, KKTMat, matArray(4)
401: ! Mat :: tmat
402: DM :: daphi, dalam
403: IS, pointer :: isglobal(:)
404: PetscErrorCode ierr
405: PetscInt its, N1, N2, i, j, irow, row(1)
406: PetscInt col(1), low, high, lamlow, lamhigh
407: PetscBool flg
408: PetscInt nloc, nloclam
409: PetscReal, parameter :: lambda_min = 0.0, lambda_max = 6.81
410: type(ex73f90tmodule_type) solver
411: PetscScalar, parameter :: one = 1.0
412: PetscScalar bval(1), cval(1)
413: PetscBool useobjective
415: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
416: ! Initialize program
417: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
418: PetscCallA(PetscInitialize(ierr))
419: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, solver%rank, ierr))
421: ! Initialize problem parameters
422: solver%lambda = 6.0
423: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', solver%lambda, flg, ierr))
424: PetscCheckA(solver%lambda <= lambda_max .and. solver%lambda >= lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')
425: useobjective = PETSC_FALSE
426: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-objective', useobjective, flg, ierr))
428: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429: ! Create vector data structures; set function evaluation routine
430: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: ! just get size
433: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 4_PETSC_INT_KIND, 4_PETSC_INT_KIND, PETSC_DECIDE, PETSC_DECIDE, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, daphi, ierr))
434: PetscCallA(DMSetFromOptions(daphi, ierr))
435: PetscCallA(DMSetUp(daphi, ierr))
436: 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))
437: N1 = solver%my*solver%mx
438: N2 = solver%my
439: flg = .false.
440: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_constraints', flg, flg, ierr))
441: if (flg) then
442: N2 = 0
443: end if
445: PetscCallA(DMDestroy(daphi, ierr))
447: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
448: ! Create matrix data structure; set Jacobian evaluation routine
449: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
450: PetscCallA(DMShellCreate(PETSC_COMM_WORLD, daphi, ierr))
451: PetscCallA(DMSetOptionsPrefix(daphi, 'phi_', ierr))
452: PetscCallA(DMSetFromOptions(daphi, ierr))
454: PetscCallA(VecCreate(PETSC_COMM_WORLD, x1, ierr))
455: PetscCallA(VecSetSizes(x1, PETSC_DECIDE, N1, ierr))
456: PetscCallA(VecSetFromOptions(x1, ierr))
458: PetscCallA(VecGetOwnershipRange(x1, low, high, ierr))
459: nloc = high - low
461: PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr))
462: PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
463: PetscCallA(MatSetUp(Amat, ierr))
465: PetscCallA(MatCreate(PETSC_COMM_WORLD, solver%AmatLin, ierr))
466: PetscCallA(MatSetSizes(solver%AmatLin, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
467: PetscCallA(MatSetUp(solver%AmatLin, ierr))
469: PetscCallA(FormJacobianLocal(x1, solver%AmatLin, solver, .false., ierr))
470: PetscCallA(MatAssemblyBegin(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
471: PetscCallA(MatAssemblyEnd(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
473: PetscCallA(DMShellSetGlobalVector(daphi, x1, ierr))
474: PetscCallA(DMShellSetMatrix(daphi, Amat, ierr))
476: PetscCallA(VecCreate(PETSC_COMM_SELF, x1loc, ierr))
477: PetscCallA(VecSetSizes(x1loc, nloc, nloc, ierr))
478: PetscCallA(VecSetFromOptions(x1loc, ierr))
479: PetscCallA(DMShellSetLocalVector(daphi, x1loc, ierr))
481: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
482: ! Create B, C, & D matrices
483: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
484: PetscCallA(MatCreate(PETSC_COMM_WORLD, Cmat, ierr))
485: PetscCallA(MatSetSizes(Cmat, PETSC_DECIDE, PETSC_DECIDE, N2, N1, ierr))
486: PetscCallA(MatSetUp(Cmat, ierr))
487: ! create data for C and B
488: PetscCallA(MatCreate(PETSC_COMM_WORLD, Bmat, ierr))
489: PetscCallA(MatSetSizes(Bmat, PETSC_DECIDE, PETSC_DECIDE, N1, N2, ierr))
490: PetscCallA(MatSetUp(Bmat, ierr))
491: ! create data for D
492: PetscCallA(MatCreate(PETSC_COMM_WORLD, Dmat, ierr))
493: PetscCallA(MatSetSizes(Dmat, PETSC_DECIDE, PETSC_DECIDE, N2, N2, ierr))
494: PetscCallA(MatSetUp(Dmat, ierr))
496: PetscCallA(VecCreate(PETSC_COMM_WORLD, x2, ierr))
497: PetscCallA(VecSetSizes(x2, PETSC_DECIDE, N2, ierr))
498: PetscCallA(VecSetFromOptions(x2, ierr))
500: PetscCallA(VecGetOwnershipRange(x2, lamlow, lamhigh, ierr))
501: nloclam = lamhigh - lamlow
503: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
504: ! Set fake B and C
505: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506: if (N2 > 0) then
507: bval(1) = -one/(solver%mx - 2)
508: ! cval = -one/(solver%my*solver%mx)
509: cval(1) = -one
510: do irow = low, high - 1
511: j = irow/solver%mx ! row in domain
512: i = mod(irow, solver%mx)
513: row(1) = irow
514: col(1) = j
515: if (.not. (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1)) then
516: PetscCallA(MatSetValues(Bmat, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, col, bval, INSERT_VALUES, ierr))
517: end if
518: row(1) = j
519: PetscCallA(MatSetValues(Cmat, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, row, cval, INSERT_VALUES, ierr))
520: end do
521: end if
522: PetscCallA(MatAssemblyBegin(Bmat, MAT_FINAL_ASSEMBLY, ierr))
523: PetscCallA(MatAssemblyEnd(Bmat, MAT_FINAL_ASSEMBLY, ierr))
524: PetscCallA(MatAssemblyBegin(Cmat, MAT_FINAL_ASSEMBLY, ierr))
525: PetscCallA(MatAssemblyEnd(Cmat, MAT_FINAL_ASSEMBLY, ierr))
527: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
528: ! Set D (identity)
529: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
530: do j = lamlow, lamhigh - 1
531: row(1) = j
532: cval(1) = one
533: PetscCallA(MatSetValues(Dmat, 1_PETSC_INT_KIND, row, 1_PETSC_INT_KIND, row, cval, INSERT_VALUES, ierr))
534: end do
535: PetscCallA(MatAssemblyBegin(Dmat, MAT_FINAL_ASSEMBLY, ierr))
536: PetscCallA(MatAssemblyEnd(Dmat, MAT_FINAL_ASSEMBLY, ierr))
538: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
539: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
540: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
541: PetscCallA(DMShellCreate(PETSC_COMM_WORLD, dalam, ierr))
542: PetscCallA(DMShellSetGlobalVector(dalam, x2, ierr))
543: PetscCallA(DMShellSetMatrix(dalam, Dmat, ierr))
545: PetscCallA(VecCreate(PETSC_COMM_SELF, x2loc, ierr))
546: PetscCallA(VecSetSizes(x2loc, nloclam, nloclam, ierr))
547: PetscCallA(VecSetFromOptions(x2loc, ierr))
548: PetscCallA(DMShellSetLocalVector(dalam, x2loc, ierr))
550: PetscCallA(DMSetOptionsPrefix(dalam, 'lambda_', ierr))
551: PetscCallA(DMSetFromOptions(dalam, ierr))
552: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
553: ! Create field split DA
554: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
555: PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD, solver%da, ierr))
556: PetscCallA(DMCompositeAddDM(solver%da, daphi, ierr))
557: PetscCallA(DMCompositeAddDM(solver%da, dalam, ierr))
558: PetscCallA(DMSetFromOptions(solver%da, ierr))
559: PetscCallA(DMSetUp(solver%da, ierr))
560: PetscCallA(DMCompositeGetGlobalISs(solver%da, isglobal, ierr))
561: solver%isPhi = isglobal(1)
562: PetscCallA(PetscObjectReference(solver%isPhi, ierr))
563: solver%isLambda = isglobal(2)
564: PetscCallA(PetscObjectReference(solver%isLambda, ierr))
566: ! cache matrices
567: solver%Amat = Amat
568: solver%Bmat = Bmat
569: solver%Cmat = Cmat
570: solver%Dmat = Dmat
572: matArray(1) = Amat
573: matArray(2) = Bmat
574: matArray(3) = Cmat
575: matArray(4) = Dmat
577: PetscCallA(MatCreateNest(PETSC_COMM_WORLD, 2_PETSC_INT_KIND, isglobal, 2_PETSC_INT_KIND, isglobal, matArray, KKTmat, ierr))
578: PetscCallA(DMCompositeRestoreGlobalISs(solver%da, isglobal, ierr))
579: PetscCallA(MatSetFromOptions(KKTmat, ierr))
581: ! Extract global and local vectors from DMDA; then duplicate for remaining
582: ! vectors that are the same types
583: PetscCallA(MatCreateVecs(KKTmat, x, PETSC_NULL_VEC, ierr))
584: PetscCallA(VecDuplicate(x, r, ierr))
586: PetscCallA(SNESCreate(PETSC_COMM_WORLD, mysnes, ierr))
588: PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
590: PetscCallA(SNESSetApplicationContext(mysnes, solver, ierr))
592: PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
594: ! Set function evaluation routine and vector
595: PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver, ierr))
596: if (useobjective .eqv. PETSC_TRUE) then
597: PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
598: end if
599: PetscCallA(SNESSetJacobian(mysnes, KKTmat, KKTmat, FormJacobian, solver, ierr))
601: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
602: ! Customize nonlinear solver; set runtime options
603: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
604: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
605: PetscCallA(SNESSetFromOptions(mysnes, ierr))
607: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
608: ! Evaluate initial guess; then solve nonlinear system.
609: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
610: ! Note: The user should initialize the vector, x, with the initial guess
611: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
612: ! to employ an initial guess of zero, the user should explicitly set
613: ! this vector to zero by calling VecSet().
615: PetscCallA(FormInitialGuess(mysnes, x, ierr))
616: PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr))
617: PetscCallA(SNESGetIterationNumber(mysnes, its, ierr))
618: if (solver%rank == 0) then
619: write (6, 100) its
620: end if
621: 100 format('Number of SNES iterations = ', i5)
623: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
624: ! Free work space. All PETSc objects should be destroyed when they
625: ! are no longer needed.
626: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
627: PetscCallA(MatDestroy(KKTmat, ierr))
628: PetscCallA(MatDestroy(Amat, ierr))
629: PetscCallA(MatDestroy(Dmat, ierr))
630: PetscCallA(MatDestroy(Bmat, ierr))
631: PetscCallA(MatDestroy(Cmat, ierr))
632: PetscCallA(MatDestroy(solver%AmatLin, ierr))
633: PetscCallA(ISDestroy(solver%isPhi, ierr))
634: PetscCallA(ISDestroy(solver%isLambda, ierr))
635: PetscCallA(VecDestroy(x, ierr))
636: PetscCallA(VecDestroy(x2, ierr))
637: PetscCallA(VecDestroy(x1, ierr))
638: PetscCallA(VecDestroy(x1loc, ierr))
639: PetscCallA(VecDestroy(x2loc, ierr))
640: PetscCallA(VecDestroy(r, ierr))
641: PetscCallA(SNESDestroy(mysnes, ierr))
642: PetscCallA(DMDestroy(solver%da, ierr))
643: PetscCallA(DMDestroy(daphi, ierr))
644: PetscCallA(DMDestroy(dalam, ierr))
646: PetscCallA(PetscFinalize(ierr))
647: end
649: !/*TEST
650: !
651: ! build:
652: ! requires: !single !complex
653: !
654: ! test:
655: ! nsize: 4
656: ! 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.
657: !
658: ! test:
659: ! args: -snes_linesearch_type {{secant cp}separate output} -objective {{false true}shared output}
660: !
661: ! test:
662: ! args: -snes_linesearch_type bt -objective {{false true}separate output}
663: !
664: !TEST*/