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/petscdm.h>
43: #include <petsc/finclude/petscmat.h>
44: use petscdm
45: use petscmat
46: type ex73f90tmodule_type
47: DM::da
48: ! temp A block stuff
49: PetscInt mx,my
50: PetscMPIInt rank
51: PetscReal lambda
52: ! Mats
53: Mat::Amat,AmatLin,Bmat,CMat,Dmat
54: IS::isPhi,isLambda
55: end type ex73f90tmodule_type
57: end module ex73f90tmodule
59: module ex73f90tmodule_interfaces
60: use ex73f90tmodule
62: Interface SNESSetApplicationContext
63: Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
64: #include <petsc/finclude/petscsnes.h>
65: use petscsnes
66: use ex73f90tmodule
67: SNES:: snesIn
68: type(ex73f90tmodule_type) ctx
69: PetscErrorCode ierr
70: End Subroutine
71: End Interface SNESSetApplicationContext
73: Interface SNESGetApplicationContext
74: Subroutine SNESGetApplicationContext(snesIn,ctx,ierr)
75: #include <petsc/finclude/petscsnes.h>
76: use petscsnes
77: use ex73f90tmodule
78: SNES:: snesIn
79: type(ex73f90tmodule_type), pointer :: ctx
80: PetscErrorCode ierr
81: End Subroutine
82: End Interface SNESGetApplicationContext
83: end module ex73f90tmodule_interfaces
85: program main
86: #include <petsc/finclude/petscdm.h>
87: #include <petsc/finclude/petscsnes.h>
88: use petscdm
89: use petscdmda
90: use petscsnes
91: use ex73f90tmodule
92: use ex73f90tmodule_interfaces
93: implicit none
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Variable declarations
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: !
98: ! Variables:
99: ! mysnes - nonlinear solver
100: ! x, r - solution, residual vectors
101: ! J - Jacobian matrix
102: ! its - iterations for convergence
103: ! Nx, Ny - number of preocessors in x- and y- directions
104: !
105: SNES:: mysnes
106: Vec:: x,r,x2,x1,x1loc,x2loc
107: Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
108: ! Mat:: tmat
109: DM:: daphi,dalam
110: IS:: isglobal(2)
111: PetscErrorCode ierr
112: PetscInt its,N1,N2,i,j,irow,row(1)
113: PetscInt col(1),low,high,lamlow,lamhigh
114: PetscBool flg
115: PetscInt ione,nfour,itwo,nloc,nloclam
116: PetscReal lambda_max,lambda_min
117: type(ex73f90tmodule_type) solver
118: PetscScalar bval(1),cval(1),one
120: ! Note: Any user-defined Fortran routines (such as FormJacobian)
121: ! MUST be declared as external.
122: external FormInitialGuess,FormJacobian,FormFunction
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Initialize program
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: PetscCallA(PetscInitialize(ierr))
128: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr))
130: ! Initialize problem parameters
131: lambda_max = 6.81_PETSC_REAL_KIND
132: lambda_min = 0.0
133: solver%lambda = 6.0
134: ione = 1
135: nfour = 4
136: itwo = 2
137: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr))
138: PetscCheckA(solver%lambda .le. lambda_max .and. solver%lambda .ge. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range')
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: ! Create vector data structures; set function evaluation routine
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! just get size
145: 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,PETSC_NULL_INTEGER,daphi,ierr))
146: PetscCallA(DMSetFromOptions(daphi,ierr))
147: PetscCallA(DMSetUp(daphi,ierr))
148: 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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
149: N1 = solver%my*solver%mx
150: N2 = solver%my
151: flg = .false.
152: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr))
153: if (flg) then
154: N2 = 0
155: endif
157: PetscCallA(DMDestroy(daphi,ierr))
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Create matrix data structure; set Jacobian evaluation routine
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,daphi,ierr))
163: PetscCallA(DMSetOptionsPrefix(daphi,'phi_',ierr))
164: PetscCallA(DMSetFromOptions(daphi,ierr))
166: PetscCallA(VecCreate(PETSC_COMM_WORLD,x1,ierr))
167: PetscCallA(VecSetSizes(x1,PETSC_DECIDE,N1,ierr))
168: PetscCallA(VecSetFromOptions(x1,ierr))
170: PetscCallA(VecGetOwnershipRange(x1,low,high,ierr))
171: nloc = high - low
173: PetscCallA(MatCreate(PETSC_COMM_WORLD,Amat,ierr))
174: PetscCallA(MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
175: PetscCallA(MatSetUp(Amat,ierr))
177: PetscCallA(MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr))
178: PetscCallA(MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
179: PetscCallA(MatSetUp(solver%AmatLin,ierr))
181: PetscCallA(FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr))
182: PetscCallA(MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
183: PetscCallA(MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
185: PetscCallA(DMShellSetGlobalVector(daphi,x1,ierr))
186: PetscCallA(DMShellSetMatrix(daphi,Amat,ierr))
188: PetscCallA(VecCreate(PETSC_COMM_SELF,x1loc,ierr))
189: PetscCallA(VecSetSizes(x1loc,nloc,nloc,ierr))
190: PetscCallA(VecSetFromOptions(x1loc,ierr))
191: PetscCallA(DMShellSetLocalVector(daphi,x1loc,ierr))
193: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: ! Create B, C, & D matrices
195: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: PetscCallA(MatCreate(PETSC_COMM_WORLD,Cmat,ierr))
197: PetscCallA(MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr))
198: PetscCallA(MatSetUp(Cmat,ierr))
199: ! create data for C and B
200: PetscCallA(MatCreate(PETSC_COMM_WORLD,Bmat,ierr))
201: PetscCallA(MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr))
202: PetscCallA(MatSetUp(Bmat,ierr))
203: ! create data for D
204: PetscCallA(MatCreate(PETSC_COMM_WORLD,Dmat,ierr))
205: PetscCallA(MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr))
206: PetscCallA(MatSetUp(Dmat,ierr))
208: PetscCallA(VecCreate(PETSC_COMM_WORLD,x2,ierr))
209: PetscCallA(VecSetSizes(x2,PETSC_DECIDE,N2,ierr))
210: PetscCallA(VecSetFromOptions(x2,ierr))
212: PetscCallA(VecGetOwnershipRange(x2,lamlow,lamhigh,ierr))
213: nloclam = lamhigh-lamlow
215: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: ! Set fake B and C
217: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: one = 1.0
219: if (N2 .gt. 0) then
220: bval(1) = -one/(solver%mx-2)
221: ! cval = -one/(solver%my*solver%mx)
222: cval(1) = -one
223: do 20 irow=low,high-1
224: j = irow/solver%mx ! row in domain
225: i = mod(irow,solver%mx)
226: row(1) = irow
227: col(1) = j
228: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
229: ! no op
230: else
231: PetscCallA(MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr))
232: endif
233: row(1) = j
234: PetscCallA(MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
235: 20 continue
236: endif
237: PetscCallA(MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr))
238: PetscCallA(MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr))
239: PetscCallA(MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr))
240: PetscCallA(MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr))
242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: ! Set D (identity)
244: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: do 30 j=lamlow,lamhigh-1
246: row(1) = j
247: cval(1) = one
248: PetscCallA(MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
249: 30 continue
250: PetscCallA(MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr))
251: PetscCallA(MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr))
253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
255: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,dalam,ierr))
257: PetscCallA(DMShellSetGlobalVector(dalam,x2,ierr))
258: PetscCallA(DMShellSetMatrix(dalam,Dmat,ierr))
260: PetscCallA(VecCreate(PETSC_COMM_SELF,x2loc,ierr))
261: PetscCallA(VecSetSizes(x2loc,nloclam,nloclam,ierr))
262: PetscCallA(VecSetFromOptions(x2loc,ierr))
263: PetscCallA(DMShellSetLocalVector(dalam,x2loc,ierr))
265: PetscCallA(DMSetOptionsPrefix(dalam,'lambda_',ierr))
266: PetscCallA(DMSetFromOptions(dalam,ierr))
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: ! Create field split DA
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr))
271: PetscCallA(DMCompositeAddDM(solver%da,daphi,ierr))
272: PetscCallA(DMCompositeAddDM(solver%da,dalam,ierr))
273: PetscCallA(DMSetFromOptions(solver%da,ierr))
274: PetscCallA(DMSetUp(solver%da,ierr))
275: PetscCallA(DMCompositeGetGlobalISs(solver%da,isglobal,ierr))
276: solver%isPhi = isglobal(1)
277: solver%isLambda = isglobal(2)
279: ! cache matrices
280: solver%Amat = Amat
281: solver%Bmat = Bmat
282: solver%Cmat = Cmat
283: solver%Dmat = Dmat
285: matArray(1) = Amat
286: matArray(2) = Bmat
287: matArray(3) = Cmat
288: matArray(4) = Dmat
290: PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr))
291: PetscCallA(MatSetFromOptions(KKTmat,ierr))
293: ! Extract global and local vectors from DMDA; then duplicate for remaining
294: ! vectors that are the same types
295: PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr))
296: PetscCallA(VecDuplicate(x,r,ierr))
298: PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr))
300: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
302: PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr))
304: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
306: ! Set function evaluation routine and vector
307: PetscCallA(SNESSetFunction(mysnes,r,FormFunction,solver,ierr))
309: PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr))
311: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
312: ! Customize nonlinear solver; set runtime options
313: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
315: PetscCallA(SNESSetFromOptions(mysnes,ierr))
317: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318: ! Evaluate initial guess; then solve nonlinear system.
319: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: ! Note: The user should initialize the vector, x, with the initial guess
321: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
322: ! to employ an initial guess of zero, the user should explicitly set
323: ! this vector to zero by calling VecSet().
325: PetscCallA(FormInitialGuess(mysnes,x,ierr))
326: PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr))
327: PetscCallA(SNESGetIterationNumber(mysnes,its,ierr))
328: if (solver%rank .eq. 0) then
329: write(6,100) its
330: endif
331: 100 format('Number of SNES iterations = ',i5)
333: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
334: ! Free work space. All PETSc objects should be destroyed when they
335: ! are no longer needed.
336: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337: PetscCallA(MatDestroy(KKTmat,ierr))
338: PetscCallA(MatDestroy(Amat,ierr))
339: PetscCallA(MatDestroy(Dmat,ierr))
340: PetscCallA(MatDestroy(Bmat,ierr))
341: PetscCallA(MatDestroy(Cmat,ierr))
342: PetscCallA(MatDestroy(solver%AmatLin,ierr))
343: PetscCallA(ISDestroy(solver%isPhi,ierr))
344: PetscCallA(ISDestroy(solver%isLambda,ierr))
345: PetscCallA(VecDestroy(x,ierr))
346: PetscCallA(VecDestroy(x2,ierr))
347: PetscCallA(VecDestroy(x1,ierr))
348: PetscCallA(VecDestroy(x1loc,ierr))
349: PetscCallA(VecDestroy(x2loc,ierr))
350: PetscCallA(VecDestroy(r,ierr))
351: PetscCallA(SNESDestroy(mysnes,ierr))
352: PetscCallA(DMDestroy(solver%da,ierr))
353: PetscCallA(DMDestroy(daphi,ierr))
354: PetscCallA(DMDestroy(dalam,ierr))
356: PetscCallA(PetscFinalize(ierr))
357: end
359: ! ---------------------------------------------------------------------
360: !
361: ! FormInitialGuess - Forms initial approximation.
362: !
363: ! Input Parameters:
364: ! X - vector
365: !
366: ! Output Parameter:
367: ! X - vector
368: !
369: ! Notes:
370: ! This routine serves as a wrapper for the lower-level routine
371: ! "InitialGuessLocal", where the actual computations are
372: ! done using the standard Fortran style of treating the local
373: ! vector data as a multidimensional array over the local mesh.
374: ! This routine merely handles ghost point scatters and accesses
375: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
376: !
377: subroutine FormInitialGuess(mysnes,Xnest,ierr)
378: #include <petsc/finclude/petscsnes.h>
379: use petscsnes
380: use ex73f90tmodule
381: use ex73f90tmodule_interfaces
382: implicit none
383: ! Input/output variables:
384: SNES:: mysnes
385: Vec:: Xnest
386: PetscErrorCode ierr
388: ! Declarations for use with local arrays:
389: type(ex73f90tmodule_type), pointer:: solver
390: Vec:: Xsub(2)
391: PetscInt:: izero,ione,itwo
393: izero = 0
394: ione = 1
395: itwo = 2
396: ierr = 0
397: PetscCall(SNESGetApplicationContext(mysnes,solver,ierr))
398: PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
400: PetscCall(InitialGuessLocal(solver,Xsub(1),ierr))
401: PetscCall(VecAssemblyBegin(Xsub(1),ierr))
402: PetscCall(VecAssemblyEnd(Xsub(1),ierr))
404: ! zero out lambda
405: PetscCall(VecZeroEntries(Xsub(2),ierr))
406: PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
408: end subroutine FormInitialGuess
410: ! ---------------------------------------------------------------------
411: !
412: ! InitialGuessLocal - Computes initial approximation, called by
413: ! the higher level routine FormInitialGuess().
414: !
415: ! Input Parameter:
416: ! X1 - local vector data
417: !
418: ! Output Parameters:
419: ! x - local vector data
420: ! ierr - error code
421: !
422: ! Notes:
423: ! This routine uses standard Fortran-style computations over a 2-dim array.
424: !
425: subroutine InitialGuessLocal(solver,X1,ierr)
426: #include <petsc/finclude/petscsys.h>
427: use petscsys
428: use ex73f90tmodule
429: implicit none
430: ! Input/output variables:
431: type (ex73f90tmodule_type) solver
432: Vec:: X1
433: PetscErrorCode ierr
435: ! Local variables:
436: PetscInt row,i,j,ione,low,high
437: PetscReal temp1,temp,hx,hy,v
438: PetscReal one
440: ! Set parameters
441: ione = 1
442: ierr = 0
443: one = 1.0
444: hx = one/(solver%mx-1)
445: hy = one/(solver%my-1)
446: temp1 = solver%lambda/(solver%lambda + one) + one
448: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
450: do 20 row=low,high-1
451: j = row/solver%mx
452: i = mod(row,solver%mx)
453: temp = min(j,solver%my-j+1)*hy
454: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
455: v = 0.0
456: else
457: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
458: endif
459: PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr))
460: 20 continue
462: end subroutine InitialGuessLocal
464: ! ---------------------------------------------------------------------
465: !
466: ! FormJacobian - Evaluates Jacobian matrix.
467: !
468: ! Input Parameters:
469: ! dummy - the SNES context
470: ! x - input vector
471: ! solver - solver data
472: !
473: ! Output Parameters:
474: ! jac - Jacobian matrix
475: ! jac_prec - optionally different preconditioning matrix (not used here)
476: ! flag - flag indicating matrix structure
477: !
478: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
479: #include <petsc/finclude/petscsnes.h>
480: use petscsnes
481: use ex73f90tmodule
482: implicit none
483: ! Input/output variables:
484: SNES:: dummy
485: Vec:: X
486: Mat:: jac,jac_prec
487: type(ex73f90tmodule_type) solver
488: PetscErrorCode ierr
490: ! Declarations for use with local arrays:
491: Vec:: Xsub(1)
492: Mat:: Amat
493: PetscInt ione
495: ione = 1
497: PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
499: ! Compute entries for the locally owned part of the Jacobian preconditioner.
500: PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr))
502: PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr))
503: PetscCall(MatDestroy(Amat,ierr)) ! discard our reference
504: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
506: ! the rest of the matrix is not touched
507: PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
508: PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
509: if (jac .ne. jac_prec) then
510: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
511: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
512: end if
514: ! Tell the matrix we will never add a new nonzero location to the
515: ! matrix. If we do it will generate an error.
516: PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
518: end subroutine FormJacobian
520: ! ---------------------------------------------------------------------
521: !
522: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
523: ! called by the higher level routine FormJacobian().
524: !
525: ! Input Parameters:
526: ! x - local vector data
527: !
528: ! Output Parameters:
529: ! jac - Jacobian preconditioner matrix
530: ! ierr - error code
531: !
532: ! Notes:
533: ! This routine uses standard Fortran-style computations over a 2-dim array.
534: !
535: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
536: #include <petsc/finclude/petscmat.h>
537: use petscmat
538: use ex73f90tmodule
539: implicit none
540: ! Input/output variables:
541: type (ex73f90tmodule_type) solver
542: Vec:: X1
543: Mat:: jac
544: logical add_nl_term
545: PetscErrorCode ierr
547: ! Local variables:
548: PetscInt irow,row(1),col(5),i,j
549: PetscInt ione,ifive,low,high,ii
550: PetscScalar two,one,hx,hy,hy2inv
551: PetscScalar hx2inv,sc,v(5)
552: PetscScalar,pointer :: lx_v(:)
554: ! Set parameters
555: ione = 1
556: ifive = 5
557: one = 1.0
558: two = 2.0
559: hx = one/(solver%mx-1)
560: hy = one/(solver%my-1)
561: sc = solver%lambda
562: hx2inv = one/(hx*hx)
563: hy2inv = one/(hy*hy)
565: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
566: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
568: ii = 0
569: do 20 irow=low,high-1
570: j = irow/solver%mx
571: i = mod(irow,solver%mx)
572: ii = ii + 1 ! one based local index
573: ! boundary points
574: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
575: col(1) = irow
576: row(1) = irow
577: v(1) = one
578: PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr))
579: ! interior grid points
580: else
581: v(1) = -hy2inv
582: if (j-1==0) v(1) = 0.0
583: v(2) = -hx2inv
584: if (i-1==0) v(2) = 0.0
585: v(3) = two*(hx2inv + hy2inv)
586: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
587: v(4) = -hx2inv
588: if (i+1==solver%mx-1) v(4) = 0.0
589: v(5) = -hy2inv
590: if (j+1==solver%my-1) v(5) = 0.0
591: col(1) = irow - solver%mx
592: col(2) = irow - 1
593: col(3) = irow
594: col(4) = irow + 1
595: col(5) = irow + solver%mx
596: row(1) = irow
597: PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr))
598: endif
599: 20 continue
601: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
603: end subroutine FormJacobianLocal
605: ! ---------------------------------------------------------------------
606: !
607: ! FormFunction - Evaluates nonlinear function, F(x).
608: !
609: ! Input Parameters:
610: ! snes - the SNES context
611: ! X - input vector
612: ! dummy - optional user-defined context, as set by SNESSetFunction()
613: ! (not used here)
614: !
615: ! Output Parameter:
616: ! F - function vector
617: !
618: subroutine FormFunction(snesIn,X,F,solver,ierr)
619: #include <petsc/finclude/petscsnes.h>
620: use petscsnes
621: use ex73f90tmodule
622: implicit none
623: ! Input/output variables:
624: SNES:: snesIn
625: Vec:: X,F
626: PetscErrorCode ierr
627: type (ex73f90tmodule_type) solver
629: ! Declarations for use with local arrays:
630: Vec:: Xsub(2),Fsub(2)
631: PetscInt itwo
633: ! Scatter ghost points to local vector, using the 2-step process
634: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
635: ! By placing code between these two statements, computations can
636: ! be done while messages are in transition.
638: itwo = 2
639: PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
640: PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
642: PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr))
643: PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
645: ! do rest of operator (linear)
646: PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr))
647: PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
648: PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
650: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
651: PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
652: end subroutine formfunction
654: ! ---------------------------------------------------------------------
655: !
656: ! FormFunctionNLTerm - Computes nonlinear function, called by
657: ! the higher level routine FormFunction().
658: !
659: ! Input Parameter:
660: ! x - local vector data
661: !
662: ! Output Parameters:
663: ! f - local vector data, f(x)
664: ! ierr - error code
665: !
666: ! Notes:
667: ! This routine uses standard Fortran-style computations over a 2-dim array.
668: !
669: subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
670: #include <petsc/finclude/petscvec.h>
671: use petscvec
672: use ex73f90tmodule
673: implicit none
674: ! Input/output variables:
675: type (ex73f90tmodule_type) solver
676: Vec:: X1,F1
677: PetscErrorCode ierr
678: ! Local variables:
679: PetscScalar sc
680: PetscScalar u,v(1)
681: PetscInt i,j,low,high,ii,ione,irow,row(1)
682: PetscScalar,pointer :: lx_v(:)
684: sc = solver%lambda
685: ione = 1
687: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
688: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
690: ! Compute function over the locally owned part of the grid
691: ii = 0
692: do 20 irow=low,high-1
693: j = irow/solver%mx
694: i = mod(irow,solver%mx)
695: ii = ii + 1 ! one based local index
696: row(1) = irow
697: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
698: v(1) = 0.0
699: else
700: u = lx_v(ii)
701: v(1) = -sc*exp(u)
702: endif
703: PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr))
704: 20 continue
706: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
708: PetscCall(VecAssemblyBegin(F1,ierr))
709: PetscCall(VecAssemblyEnd(F1,ierr))
711: ierr = 0
712: end subroutine FormFunctionNLTerm
714: !/*TEST
715: !
716: ! build:
717: ! requires: !single !complex
718: !
719: ! test:
720: ! nsize: 4
721: ! 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.
722: !
723: !TEST*/