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: if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then
139: SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range')
140: endif
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: ! Create vector data structures; set function evaluation routine
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: ! just get size
147: 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))
148: PetscCallA(DMSetFromOptions(daphi,ierr))
149: PetscCallA(DMSetUp(daphi,ierr))
150: 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))
151: N1 = solver%my*solver%mx
152: N2 = solver%my
153: flg = .false.
154: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr))
155: if (flg) then
156: N2 = 0
157: endif
159: PetscCallA(DMDestroy(daphi,ierr))
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! Create matrix data structure; set Jacobian evaluation routine
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,daphi,ierr))
165: PetscCallA(DMSetOptionsPrefix(daphi,'phi_',ierr))
166: PetscCallA(DMSetFromOptions(daphi,ierr))
168: PetscCallA(VecCreate(PETSC_COMM_WORLD,x1,ierr))
169: PetscCallA(VecSetSizes(x1,PETSC_DECIDE,N1,ierr))
170: PetscCallA(VecSetFromOptions(x1,ierr))
172: PetscCallA(VecGetOwnershipRange(x1,low,high,ierr))
173: nloc = high - low
175: PetscCallA(MatCreate(PETSC_COMM_WORLD,Amat,ierr))
176: PetscCallA(MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
177: PetscCallA(MatSetUp(Amat,ierr))
179: PetscCallA(MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr))
180: PetscCallA(MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
181: PetscCallA(MatSetUp(solver%AmatLin,ierr))
183: PetscCallA(FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr))
184: PetscCallA(MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
185: PetscCallA(MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
187: PetscCallA(DMShellSetGlobalVector(daphi,x1,ierr))
188: PetscCallA(DMShellSetMatrix(daphi,Amat,ierr))
190: PetscCallA(VecCreate(PETSC_COMM_SELF,x1loc,ierr))
191: PetscCallA(VecSetSizes(x1loc,nloc,nloc,ierr))
192: PetscCallA(VecSetFromOptions(x1loc,ierr))
193: PetscCallA(DMShellSetLocalVector(daphi,x1loc,ierr))
195: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: ! Create B, C, & D matrices
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: PetscCallA(MatCreate(PETSC_COMM_WORLD,Cmat,ierr))
199: PetscCallA(MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr))
200: PetscCallA(MatSetUp(Cmat,ierr))
201: ! create data for C and B
202: PetscCallA(MatCreate(PETSC_COMM_WORLD,Bmat,ierr))
203: PetscCallA(MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr))
204: PetscCallA(MatSetUp(Bmat,ierr))
205: ! create data for D
206: PetscCallA(MatCreate(PETSC_COMM_WORLD,Dmat,ierr))
207: PetscCallA(MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr))
208: PetscCallA(MatSetUp(Dmat,ierr))
210: PetscCallA(VecCreate(PETSC_COMM_WORLD,x2,ierr))
211: PetscCallA(VecSetSizes(x2,PETSC_DECIDE,N2,ierr))
212: PetscCallA(VecSetFromOptions(x2,ierr))
214: PetscCallA(VecGetOwnershipRange(x2,lamlow,lamhigh,ierr))
215: nloclam = lamhigh-lamlow
217: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: ! Set fake B and C
219: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: one = 1.0
221: if (N2 .gt. 0) then
222: bval(1) = -one/(solver%mx-2)
223: ! cval = -one/(solver%my*solver%mx)
224: cval(1) = -one
225: do 20 irow=low,high-1
226: j = irow/solver%mx ! row in domain
227: i = mod(irow,solver%mx)
228: row(1) = irow
229: col(1) = j
230: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
231: ! no op
232: else
233: PetscCallA(MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr))
234: endif
235: row(1) = j
236: PetscCallA(MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
237: 20 continue
238: endif
239: PetscCallA(MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr))
240: PetscCallA(MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr))
241: PetscCallA(MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr))
242: PetscCallA(MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr))
244: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: ! Set D (identity)
246: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247: do 30 j=lamlow,lamhigh-1
248: row(1) = j
249: cval(1) = one
250: PetscCallA(MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
251: 30 continue
252: PetscCallA(MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr))
253: PetscCallA(MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr))
255: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,dalam,ierr))
259: PetscCallA(DMShellSetGlobalVector(dalam,x2,ierr))
260: PetscCallA(DMShellSetMatrix(dalam,Dmat,ierr))
262: PetscCallA(VecCreate(PETSC_COMM_SELF,x2loc,ierr))
263: PetscCallA(VecSetSizes(x2loc,nloclam,nloclam,ierr))
264: PetscCallA(VecSetFromOptions(x2loc,ierr))
265: PetscCallA(DMShellSetLocalVector(dalam,x2loc,ierr))
267: PetscCallA(DMSetOptionsPrefix(dalam,'lambda_',ierr))
268: PetscCallA(DMSetFromOptions(dalam,ierr))
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: ! Create field split DA
271: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272: PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr))
273: PetscCallA(DMCompositeAddDM(solver%da,daphi,ierr))
274: PetscCallA(DMCompositeAddDM(solver%da,dalam,ierr))
275: PetscCallA(DMSetFromOptions(solver%da,ierr))
276: PetscCallA(DMSetUp(solver%da,ierr))
277: PetscCallA(DMCompositeGetGlobalISs(solver%da,isglobal,ierr))
278: solver%isPhi = isglobal(1)
279: solver%isLambda = isglobal(2)
281: ! cache matrices
282: solver%Amat = Amat
283: solver%Bmat = Bmat
284: solver%Cmat = Cmat
285: solver%Dmat = Dmat
287: matArray(1) = Amat
288: matArray(2) = Bmat
289: matArray(3) = Cmat
290: matArray(4) = Dmat
292: PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr))
293: PetscCallA(MatSetFromOptions(KKTmat,ierr))
295: ! Extract global and local vectors from DMDA; then duplicate for remaining
296: ! vectors that are the same types
297: PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr))
298: PetscCallA(VecDuplicate(x,r,ierr))
300: PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr))
302: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
304: PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr))
306: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
308: ! Set function evaluation routine and vector
309: PetscCallA(SNESSetFunction(mysnes,r,FormFunction,solver,ierr))
311: PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr))
313: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314: ! Customize nonlinear solver; set runtime options
315: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
316: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
317: PetscCallA(SNESSetFromOptions(mysnes,ierr))
319: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: ! Evaluate initial guess; then solve nonlinear system.
321: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
322: ! Note: The user should initialize the vector, x, with the initial guess
323: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
324: ! to employ an initial guess of zero, the user should explicitly set
325: ! this vector to zero by calling VecSet().
327: PetscCallA(FormInitialGuess(mysnes,x,ierr))
328: PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr))
329: PetscCallA(SNESGetIterationNumber(mysnes,its,ierr))
330: if (solver%rank .eq. 0) then
331: write(6,100) its
332: endif
333: 100 format('Number of SNES iterations = ',i5)
335: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336: ! Free work space. All PETSc objects should be destroyed when they
337: ! are no longer needed.
338: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339: PetscCallA(MatDestroy(KKTmat,ierr))
340: PetscCallA(MatDestroy(Amat,ierr))
341: PetscCallA(MatDestroy(Dmat,ierr))
342: PetscCallA(MatDestroy(Bmat,ierr))
343: PetscCallA(MatDestroy(Cmat,ierr))
344: PetscCallA(MatDestroy(solver%AmatLin,ierr))
345: PetscCallA(ISDestroy(solver%isPhi,ierr))
346: PetscCallA(ISDestroy(solver%isLambda,ierr))
347: PetscCallA(VecDestroy(x,ierr))
348: PetscCallA(VecDestroy(x2,ierr))
349: PetscCallA(VecDestroy(x1,ierr))
350: PetscCallA(VecDestroy(x1loc,ierr))
351: PetscCallA(VecDestroy(x2loc,ierr))
352: PetscCallA(VecDestroy(r,ierr))
353: PetscCallA(SNESDestroy(mysnes,ierr))
354: PetscCallA(DMDestroy(solver%da,ierr))
355: PetscCallA(DMDestroy(daphi,ierr))
356: PetscCallA(DMDestroy(dalam,ierr))
358: PetscCallA(PetscFinalize(ierr))
359: end
361: ! ---------------------------------------------------------------------
362: !
363: ! FormInitialGuess - Forms initial approximation.
364: !
365: ! Input Parameters:
366: ! X - vector
367: !
368: ! Output Parameter:
369: ! X - vector
370: !
371: ! Notes:
372: ! This routine serves as a wrapper for the lower-level routine
373: ! "InitialGuessLocal", where the actual computations are
374: ! done using the standard Fortran style of treating the local
375: ! vector data as a multidimensional array over the local mesh.
376: ! This routine merely handles ghost point scatters and accesses
377: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
378: !
379: subroutine FormInitialGuess(mysnes,Xnest,ierr)
380: #include <petsc/finclude/petscsnes.h>
381: use petscsnes
382: use ex73f90tmodule
383: use ex73f90tmodule_interfaces
384: implicit none
385: ! Input/output variables:
386: SNES:: mysnes
387: Vec:: Xnest
388: PetscErrorCode ierr
390: ! Declarations for use with local arrays:
391: type(ex73f90tmodule_type), pointer:: solver
392: Vec:: Xsub(2)
393: PetscInt:: izero,ione,itwo
395: izero = 0
396: ione = 1
397: itwo = 2
398: 0
399: PetscCall(SNESGetApplicationContext(mysnes,solver,ierr))
400: PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
402: PetscCall(InitialGuessLocal(solver,Xsub(1),ierr))
403: PetscCall(VecAssemblyBegin(Xsub(1),ierr))
404: PetscCall(VecAssemblyEnd(Xsub(1),ierr))
406: ! zero out lambda
407: PetscCall(VecZeroEntries(Xsub(2),ierr))
408: PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
410: return
411: end subroutine FormInitialGuess
413: ! ---------------------------------------------------------------------
414: !
415: ! InitialGuessLocal - Computes initial approximation, called by
416: ! the higher level routine FormInitialGuess().
417: !
418: ! Input Parameter:
419: ! X1 - local vector data
420: !
421: ! Output Parameters:
422: ! x - local vector data
423: ! ierr - error code
424: !
425: ! Notes:
426: ! This routine uses standard Fortran-style computations over a 2-dim array.
427: !
428: subroutine InitialGuessLocal(solver,X1,ierr)
429: #include <petsc/finclude/petscsys.h>
430: use petscsys
431: use ex73f90tmodule
432: implicit none
433: ! Input/output variables:
434: type (ex73f90tmodule_type) solver
435: Vec:: X1
436: PetscErrorCode ierr
438: ! Local variables:
439: PetscInt row,i,j,ione,low,high
440: PetscReal temp1,temp,hx,hy,v
441: PetscReal one
443: ! Set parameters
444: ione = 1
445: 0
446: one = 1.0
447: hx = one/(solver%mx-1)
448: hy = one/(solver%my-1)
449: temp1 = solver%lambda/(solver%lambda + one) + one
451: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
453: do 20 row=low,high-1
454: j = row/solver%mx
455: i = mod(row,solver%mx)
456: temp = min(j,solver%my-j+1)*hy
457: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
458: v = 0.0
459: else
460: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
461: endif
462: PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr))
463: 20 continue
465: return
466: end subroutine InitialGuessLocal
468: ! ---------------------------------------------------------------------
469: !
470: ! FormJacobian - Evaluates Jacobian matrix.
471: !
472: ! Input Parameters:
473: ! dummy - the SNES context
474: ! x - input vector
475: ! solver - solver data
476: !
477: ! Output Parameters:
478: ! jac - Jacobian matrix
479: ! jac_prec - optionally different preconditioning matrix (not used here)
480: ! flag - flag indicating matrix structure
481: !
482: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
483: #include <petsc/finclude/petscsnes.h>
484: use petscsnes
485: use ex73f90tmodule
486: implicit none
487: ! Input/output variables:
488: SNES:: dummy
489: Vec:: X
490: Mat:: jac,jac_prec
491: type(ex73f90tmodule_type) solver
492: PetscErrorCode ierr
494: ! Declarations for use with local arrays:
495: Vec:: Xsub(1)
496: Mat:: Amat
497: PetscInt ione
499: ione = 1
501: PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
503: ! Compute entries for the locally owned part of the Jacobian preconditioner.
504: PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr))
506: PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr))
507: PetscCall(MatDestroy(Amat,ierr)) ! discard our reference
508: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
510: ! the rest of the matrix is not touched
511: PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
512: PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
513: if (jac .ne. jac_prec) then
514: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
515: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
516: end if
518: ! Tell the matrix we will never add a new nonzero location to the
519: ! matrix. If we do it will generate an error.
520: PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
522: return
523: end subroutine FormJacobian
525: ! ---------------------------------------------------------------------
526: !
527: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
528: ! called by the higher level routine FormJacobian().
529: !
530: ! Input Parameters:
531: ! x - local vector data
532: !
533: ! Output Parameters:
534: ! jac - Jacobian preconditioner matrix
535: ! ierr - error code
536: !
537: ! Notes:
538: ! This routine uses standard Fortran-style computations over a 2-dim array.
539: !
540: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
541: #include <petsc/finclude/petscmat.h>
542: use petscmat
543: use ex73f90tmodule
544: implicit none
545: ! Input/output variables:
546: type (ex73f90tmodule_type) solver
547: Vec:: X1
548: Mat:: jac
549: logical add_nl_term
550: PetscErrorCode ierr
552: ! Local variables:
553: PetscInt irow,row(1),col(5),i,j
554: PetscInt ione,ifive,low,high,ii
555: PetscScalar two,one,hx,hy,hy2inv
556: PetscScalar hx2inv,sc,v(5)
557: PetscScalar,pointer :: lx_v(:)
559: ! Set parameters
560: ione = 1
561: ifive = 5
562: one = 1.0
563: two = 2.0
564: hx = one/(solver%mx-1)
565: hy = one/(solver%my-1)
566: sc = solver%lambda
567: hx2inv = one/(hx*hx)
568: hy2inv = one/(hy*hy)
570: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
571: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
573: ii = 0
574: do 20 irow=low,high-1
575: j = irow/solver%mx
576: i = mod(irow,solver%mx)
577: ii = ii + 1 ! one based local index
578: ! boundary points
579: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
580: col(1) = irow
581: row(1) = irow
582: v(1) = one
583: PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr))
584: ! interior grid points
585: else
586: v(1) = -hy2inv
587: if (j-1==0) v(1) = 0.0
588: v(2) = -hx2inv
589: if (i-1==0) v(2) = 0.0
590: v(3) = two*(hx2inv + hy2inv)
591: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
592: v(4) = -hx2inv
593: if (i+1==solver%mx-1) v(4) = 0.0
594: v(5) = -hy2inv
595: if (j+1==solver%my-1) v(5) = 0.0
596: col(1) = irow - solver%mx
597: col(2) = irow - 1
598: col(3) = irow
599: col(4) = irow + 1
600: col(5) = irow + solver%mx
601: row(1) = irow
602: PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr))
603: endif
604: 20 continue
606: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
608: return
609: end subroutine FormJacobianLocal
611: ! ---------------------------------------------------------------------
612: !
613: ! FormFunction - Evaluates nonlinear function, F(x).
614: !
615: ! Input Parameters:
616: ! snes - the SNES context
617: ! X - input vector
618: ! dummy - optional user-defined context, as set by SNESSetFunction()
619: ! (not used here)
620: !
621: ! Output Parameter:
622: ! F - function vector
623: !
624: subroutine FormFunction(snesIn,X,F,solver,ierr)
625: #include <petsc/finclude/petscsnes.h>
626: use petscsnes
627: use ex73f90tmodule
628: implicit none
629: ! Input/output variables:
630: SNES:: snesIn
631: Vec:: X,F
632: PetscErrorCode ierr
633: type (ex73f90tmodule_type) solver
635: ! Declarations for use with local arrays:
636: Vec:: Xsub(2),Fsub(2)
637: PetscInt itwo
639: ! Scatter ghost points to local vector, using the 2-step process
640: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
641: ! By placing code between these two statements, computations can
642: ! be done while messages are in transition.
644: itwo = 2
645: PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
646: PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
648: PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr))
649: PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
651: ! do rest of operator (linear)
652: PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr))
653: PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
654: PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
656: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
657: PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
658: return
659: end subroutine formfunction
661: ! ---------------------------------------------------------------------
662: !
663: ! FormFunctionNLTerm - Computes nonlinear function, called by
664: ! the higher level routine FormFunction().
665: !
666: ! Input Parameter:
667: ! x - local vector data
668: !
669: ! Output Parameters:
670: ! f - local vector data, f(x)
671: ! ierr - error code
672: !
673: ! Notes:
674: ! This routine uses standard Fortran-style computations over a 2-dim array.
675: !
676: subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
677: #include <petsc/finclude/petscvec.h>
678: use petscvec
679: use ex73f90tmodule
680: implicit none
681: ! Input/output variables:
682: type (ex73f90tmodule_type) solver
683: Vec:: X1,F1
684: PetscErrorCode ierr
685: ! Local variables:
686: PetscScalar sc
687: PetscScalar u,v(1)
688: PetscInt i,j,low,high,ii,ione,irow,row(1)
689: PetscScalar,pointer :: lx_v(:)
691: sc = solver%lambda
692: ione = 1
694: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
695: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
697: ! Compute function over the locally owned part of the grid
698: ii = 0
699: do 20 irow=low,high-1
700: j = irow/solver%mx
701: i = mod(irow,solver%mx)
702: ii = ii + 1 ! one based local index
703: row(1) = irow
704: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
705: v(1) = 0.0
706: else
707: u = lx_v(ii)
708: v(1) = -sc*exp(u)
709: endif
710: PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr))
711: 20 continue
713: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
715: PetscCall(VecAssemblyBegin(F1,ierr))
716: PetscCall(VecAssemblyEnd(F1,ierr))
718: 0
719: return
720: end subroutine FormFunctionNLTerm
722: !/*TEST
723: !
724: ! build:
725: ! requires: !single !complex
726: !
727: ! test:
728: ! nsize: 4
729: ! 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.
730: !
731: !TEST*/