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: return
409: end subroutine FormInitialGuess
411: ! ---------------------------------------------------------------------
412: !
413: ! InitialGuessLocal - Computes initial approximation, called by
414: ! the higher level routine FormInitialGuess().
415: !
416: ! Input Parameter:
417: ! X1 - local vector data
418: !
419: ! Output Parameters:
420: ! x - local vector data
421: ! ierr - error code
422: !
423: ! Notes:
424: ! This routine uses standard Fortran-style computations over a 2-dim array.
425: !
426: subroutine InitialGuessLocal(solver,X1,ierr)
427: #include <petsc/finclude/petscsys.h>
428: use petscsys
429: use ex73f90tmodule
430: implicit none
431: ! Input/output variables:
432: type (ex73f90tmodule_type) solver
433: Vec:: X1
434: PetscErrorCode ierr
436: ! Local variables:
437: PetscInt row,i,j,ione,low,high
438: PetscReal temp1,temp,hx,hy,v
439: PetscReal one
441: ! Set parameters
442: ione = 1
443: ierr = 0
444: one = 1.0
445: hx = one/(solver%mx-1)
446: hy = one/(solver%my-1)
447: temp1 = solver%lambda/(solver%lambda + one) + one
449: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
451: do 20 row=low,high-1
452: j = row/solver%mx
453: i = mod(row,solver%mx)
454: temp = min(j,solver%my-j+1)*hy
455: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
456: v = 0.0
457: else
458: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
459: endif
460: PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr))
461: 20 continue
463: return
464: end subroutine InitialGuessLocal
466: ! ---------------------------------------------------------------------
467: !
468: ! FormJacobian - Evaluates Jacobian matrix.
469: !
470: ! Input Parameters:
471: ! dummy - the SNES context
472: ! x - input vector
473: ! solver - solver data
474: !
475: ! Output Parameters:
476: ! jac - Jacobian matrix
477: ! jac_prec - optionally different preconditioning matrix (not used here)
478: ! flag - flag indicating matrix structure
479: !
480: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
481: #include <petsc/finclude/petscsnes.h>
482: use petscsnes
483: use ex73f90tmodule
484: implicit none
485: ! Input/output variables:
486: SNES:: dummy
487: Vec:: X
488: Mat:: jac,jac_prec
489: type(ex73f90tmodule_type) solver
490: PetscErrorCode ierr
492: ! Declarations for use with local arrays:
493: Vec:: Xsub(1)
494: Mat:: Amat
495: PetscInt ione
497: ione = 1
499: PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
501: ! Compute entries for the locally owned part of the Jacobian preconditioner.
502: PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr))
504: PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr))
505: PetscCall(MatDestroy(Amat,ierr)) ! discard our reference
506: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr))
508: ! the rest of the matrix is not touched
509: PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
510: PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
511: if (jac .ne. jac_prec) then
512: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
513: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
514: end if
516: ! Tell the matrix we will never add a new nonzero location to the
517: ! matrix. If we do it will generate an error.
518: PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
520: return
521: end subroutine FormJacobian
523: ! ---------------------------------------------------------------------
524: !
525: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
526: ! called by the higher level routine FormJacobian().
527: !
528: ! Input Parameters:
529: ! x - local vector data
530: !
531: ! Output Parameters:
532: ! jac - Jacobian preconditioner matrix
533: ! ierr - error code
534: !
535: ! Notes:
536: ! This routine uses standard Fortran-style computations over a 2-dim array.
537: !
538: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
539: #include <petsc/finclude/petscmat.h>
540: use petscmat
541: use ex73f90tmodule
542: implicit none
543: ! Input/output variables:
544: type (ex73f90tmodule_type) solver
545: Vec:: X1
546: Mat:: jac
547: logical add_nl_term
548: PetscErrorCode ierr
550: ! Local variables:
551: PetscInt irow,row(1),col(5),i,j
552: PetscInt ione,ifive,low,high,ii
553: PetscScalar two,one,hx,hy,hy2inv
554: PetscScalar hx2inv,sc,v(5)
555: PetscScalar,pointer :: lx_v(:)
557: ! Set parameters
558: ione = 1
559: ifive = 5
560: one = 1.0
561: two = 2.0
562: hx = one/(solver%mx-1)
563: hy = one/(solver%my-1)
564: sc = solver%lambda
565: hx2inv = one/(hx*hx)
566: hy2inv = one/(hy*hy)
568: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
569: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
571: ii = 0
572: do 20 irow=low,high-1
573: j = irow/solver%mx
574: i = mod(irow,solver%mx)
575: ii = ii + 1 ! one based local index
576: ! boundary points
577: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
578: col(1) = irow
579: row(1) = irow
580: v(1) = one
581: PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr))
582: ! interior grid points
583: else
584: v(1) = -hy2inv
585: if (j-1==0) v(1) = 0.0
586: v(2) = -hx2inv
587: if (i-1==0) v(2) = 0.0
588: v(3) = two*(hx2inv + hy2inv)
589: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
590: v(4) = -hx2inv
591: if (i+1==solver%mx-1) v(4) = 0.0
592: v(5) = -hy2inv
593: if (j+1==solver%my-1) v(5) = 0.0
594: col(1) = irow - solver%mx
595: col(2) = irow - 1
596: col(3) = irow
597: col(4) = irow + 1
598: col(5) = irow + solver%mx
599: row(1) = irow
600: PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr))
601: endif
602: 20 continue
604: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
606: return
607: end subroutine FormJacobianLocal
609: ! ---------------------------------------------------------------------
610: !
611: ! FormFunction - Evaluates nonlinear function, F(x).
612: !
613: ! Input Parameters:
614: ! snes - the SNES context
615: ! X - input vector
616: ! dummy - optional user-defined context, as set by SNESSetFunction()
617: ! (not used here)
618: !
619: ! Output Parameter:
620: ! F - function vector
621: !
622: subroutine FormFunction(snesIn,X,F,solver,ierr)
623: #include <petsc/finclude/petscsnes.h>
624: use petscsnes
625: use ex73f90tmodule
626: implicit none
627: ! Input/output variables:
628: SNES:: snesIn
629: Vec:: X,F
630: PetscErrorCode ierr
631: type (ex73f90tmodule_type) solver
633: ! Declarations for use with local arrays:
634: Vec:: Xsub(2),Fsub(2)
635: PetscInt itwo
637: ! Scatter ghost points to local vector, using the 2-step process
638: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
639: ! By placing code between these two statements, computations can
640: ! be done while messages are in transition.
642: itwo = 2
643: PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
644: PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
646: PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr))
647: PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
649: ! do rest of operator (linear)
650: PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr))
651: PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
652: PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
654: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr))
655: PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr))
656: return
657: end subroutine formfunction
659: ! ---------------------------------------------------------------------
660: !
661: ! FormFunctionNLTerm - Computes nonlinear function, called by
662: ! the higher level routine FormFunction().
663: !
664: ! Input Parameter:
665: ! x - local vector data
666: !
667: ! Output Parameters:
668: ! f - local vector data, f(x)
669: ! ierr - error code
670: !
671: ! Notes:
672: ! This routine uses standard Fortran-style computations over a 2-dim array.
673: !
674: subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
675: #include <petsc/finclude/petscvec.h>
676: use petscvec
677: use ex73f90tmodule
678: implicit none
679: ! Input/output variables:
680: type (ex73f90tmodule_type) solver
681: Vec:: X1,F1
682: PetscErrorCode ierr
683: ! Local variables:
684: PetscScalar sc
685: PetscScalar u,v(1)
686: PetscInt i,j,low,high,ii,ione,irow,row(1)
687: PetscScalar,pointer :: lx_v(:)
689: sc = solver%lambda
690: ione = 1
692: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
693: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
695: ! Compute function over the locally owned part of the grid
696: ii = 0
697: do 20 irow=low,high-1
698: j = irow/solver%mx
699: i = mod(irow,solver%mx)
700: ii = ii + 1 ! one based local index
701: row(1) = irow
702: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
703: v(1) = 0.0
704: else
705: u = lx_v(ii)
706: v(1) = -sc*exp(u)
707: endif
708: PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr))
709: 20 continue
711: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
713: PetscCall(VecAssemblyBegin(F1,ierr))
714: PetscCall(VecAssemblyEnd(F1,ierr))
716: ierr = 0
717: return
718: end subroutine FormFunctionNLTerm
720: !/*TEST
721: !
722: ! build:
723: ! requires: !single !complex
724: !
725: ! test:
726: ! nsize: 4
727: ! 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.
728: !
729: !TEST*/