Actual source code: ex1f.F90
1: !
2: ! Description: This example solves a nonlinear system on 1 processor with SNES.
3: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
4: ! domain. The command line options include:
5: ! -par <parameter>, where <parameter> indicates the nonlinearity of the problem
6: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
7: ! -mx <xg>, where <xg> = number of grid points in the x-direction
8: ! -my <yg>, where <yg> = number of grid points in the y-direction
9: !
11: !
12: ! --------------------------------------------------------------------------
13: !
14: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
15: ! the partial differential equation
16: !
17: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
18: !
19: ! with boundary conditions
20: !
21: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
22: !
23: ! A finite difference approximation with the usual 5-point stencil
24: ! is used to discretize the boundary value problem to obtain a nonlinear
25: ! system of equations.
26: !
27: ! The parallel version of this code is snes/tutorials/ex5f.F
28: !
29: ! --------------------------------------------------------------------------
30: subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr)
31: #include <petsc/finclude/petscsnes.h>
32: use petscsnes
33: implicit none
34: SNES snes
35: PetscReal norm
36: Vec tmp,x,y,w
37: PetscBool changed_w,changed_y
38: PetscErrorCode ierr
39: PetscInt ctx
40: PetscScalar mone
41: MPI_Comm comm
43: character(len=PETSC_MAX_PATH_LEN) :: outputString
45: PetscCallA(PetscObjectGetComm(snes,comm,ierr))
46: PetscCallA(VecDuplicate(x,tmp,ierr))
47: mone = -1.0
48: PetscCallA(VecWAXPY(tmp,mone,x,w,ierr))
49: PetscCallA(VecNorm(tmp,NORM_2,norm,ierr))
50: PetscCallA(VecDestroy(tmp,ierr))
51: write(outputString,*) norm
52: PetscCallA(PetscPrintf(comm,'Norm of search step '//trim(outputString)//'\n',ierr))
53: end
55: program main
56: #include <petsc/finclude/petscdraw.h>
57: use petscsnes
58: implicit none
59: !
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Variable declarations
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: !
64: ! Variables:
65: ! snes - nonlinear solver
66: ! x, r - solution, residual vectors
67: ! J - Jacobian matrix
68: ! its - iterations for convergence
69: ! matrix_free - flag - 1 indicates matrix-free version
70: ! lambda - nonlinearity parameter
71: ! draw - drawing context
72: !
73: SNES snes
74: MatColoring mc
75: Vec x,r
76: PetscDraw draw
77: Mat J
78: PetscBool matrix_free,flg,fd_coloring
79: PetscErrorCode ierr
80: PetscInt its,N, mx,my,i5
81: PetscMPIInt size,rank
82: PetscReal lambda_max,lambda_min,lambda
83: MatFDColoring fdcoloring
84: ISColoring iscoloring
85: PetscBool pc
86: external postcheck
88: character(len=PETSC_MAX_PATH_LEN) :: outputString
90: PetscScalar,pointer :: lx_v(:)
92: ! Store parameters in common block
94: common /params/ lambda,mx,my,fd_coloring
96: ! Note: Any user-defined Fortran routines (such as FormJacobian)
97: ! MUST be declared as external.
99: external FormFunction,FormInitialGuess,FormJacobian
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: ! Initialize program
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: PetscCallA(PetscInitialize(ierr))
106: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
107: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
109: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')
111: ! Initialize problem parameters
112: i5 = 5
113: lambda_max = 6.81
114: lambda_min = 0.0
115: lambda = 6.0
116: mx = 4
117: my = 4
118: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
119: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
120: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
121: PetscCheckA(lambda .lt. lambda_max .and. lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range ')
122: N = mx*my
123: pc = PETSC_FALSE
124: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! Create nonlinear solver context
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
132: if (pc .eqv. PETSC_TRUE) then
133: PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
134: PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
135: endif
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: ! Create vector data structures; set function evaluation routine
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
142: PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
143: PetscCallA(VecSetFromOptions(x,ierr))
144: PetscCallA(VecDuplicate(x,r,ierr))
146: ! Set function evaluation routine and vector. Whenever the nonlinear
147: ! solver needs to evaluate the nonlinear function, it will call this
148: ! routine.
149: ! - Note that the final routine argument is the user-defined
150: ! context that provides application-specific data for the
151: ! function evaluation routine.
153: PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: ! Create matrix data structure; set Jacobian evaluation routine
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Create matrix. Here we only approximately preallocate storage space
160: ! for the Jacobian. See the users manual for a discussion of better
161: ! techniques for preallocating matrix memory.
163: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
164: if (.not. matrix_free) then
165: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER_ARRAY,J,ierr))
166: endif
168: !
169: ! This option will cause the Jacobian to be computed via finite differences
170: ! efficiently using a coloring of the columns of the matrix.
171: !
172: fd_coloring = .false.
173: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
174: if (fd_coloring) then
176: !
177: ! This initializes the nonzero structure of the Jacobian. This is artificial
178: ! because clearly if we had a routine to compute the Jacobian we won't need
179: ! to use finite differences.
180: !
181: PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
182: !
183: ! Color the matrix, i.e. determine groups of columns that share no common
184: ! rows. These columns in the Jacobian can all be computed simultaneously.
185: !
186: PetscCallA(MatColoringCreate(J,mc,ierr))
187: PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
188: PetscCallA(MatColoringSetFromOptions(mc,ierr))
189: PetscCallA(MatColoringApply(mc,iscoloring,ierr))
190: PetscCallA(MatColoringDestroy(mc,ierr))
191: !
192: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
193: ! to compute the actual Jacobians via finite differences.
194: !
195: PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
196: PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
197: PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
198: PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
199: !
200: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
201: ! to compute Jacobians.
202: !
203: PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
204: PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr))
205: PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,fdcoloring,ierr))
206: PetscCallA(ISColoringDestroy(iscoloring,ierr))
208: else if (.not. matrix_free) then
210: ! Set Jacobian matrix data structure and default Jacobian evaluation
211: ! routine. Whenever the nonlinear solver needs to compute the
212: ! Jacobian matrix, it will call this routine.
213: ! - Note that the final routine argument is the user-defined
214: ! context that provides application-specific data for the
215: ! Jacobian evaluation routine.
216: ! - The user can override with:
217: ! -snes_fd : default finite differencing approximation of Jacobian
218: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
219: ! (unless user explicitly sets preconditioner)
220: ! -snes_mf_operator : form preconditioning matrix as set by the user,
221: ! but use matrix-free approx for Jacobian-vector
222: ! products within Newton-Krylov method
223: !
224: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
225: endif
227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: ! Customize nonlinear solver; set runtime options
229: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
233: PetscCallA(SNESSetFromOptions(snes,ierr))
235: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: ! Evaluate initial guess; then solve nonlinear system.
237: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: ! Note: The user should initialize the vector, x, with the initial guess
240: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
241: ! to employ an initial guess of zero, the user should explicitly set
242: ! this vector to zero by calling VecSet().
244: PetscCallA(FormInitialGuess(x,ierr))
245: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
246: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
247: write(outputString,*) its
248: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,'Number of SNES iterations = '//trim(outputString)//'\n',ierr))
250: ! PetscDraw contour plot of solution
252: PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
253: PetscCallA(PetscDrawSetFromOptions(draw,ierr))
255: PetscCallA(VecGetArrayReadF90(x,lx_v,ierr))
256: PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v,ierr))
257: PetscCallA(VecRestoreArrayReadF90(x,lx_v,ierr))
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: ! Free work space. All PETSc objects should be destroyed when they
261: ! are no longer needed.
262: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
265: if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))
267: PetscCallA(VecDestroy(x,ierr))
268: PetscCallA(VecDestroy(r,ierr))
269: PetscCallA(SNESDestroy(snes,ierr))
270: PetscCallA(PetscDrawDestroy(draw,ierr))
271: PetscCallA(PetscFinalize(ierr))
272: end
274: ! ---------------------------------------------------------------------
275: !
276: ! FormInitialGuess - Forms initial approximation.
277: !
278: ! Input Parameter:
279: ! X - vector
280: !
281: ! Output Parameters:
282: ! X - vector
283: ! ierr - error code
284: !
285: ! Notes:
286: ! This routine serves as a wrapper for the lower-level routine
287: ! "ApplicationInitialGuess", where the actual computations are
288: ! done using the standard Fortran style of treating the local
289: ! vector data as a multidimensional array over the local mesh.
290: ! This routine merely accesses the local vector data via
291: ! VecGetArrayF90() and VecRestoreArrayF90().
292: !
293: subroutine FormInitialGuess(X,ierr)
294: use petscsnes
295: implicit none
297: ! Input/output variables:
298: Vec X
299: PetscErrorCode ierr
301: ! Declarations for use with local arrays:
302: PetscScalar,pointer :: lx_v(:)
304: ierr = 0
306: ! Get a pointer to vector data.
307: ! - VecGetArrayF90() returns a pointer to the data array.
308: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
309: ! the array.
311: PetscCallA(VecGetArrayF90(X,lx_v,ierr))
313: ! Compute initial guess
315: PetscCallA(ApplicationInitialGuess(lx_v,ierr))
317: ! Restore vector
319: PetscCallA(VecRestoreArrayF90(X,lx_v,ierr))
321: end
323: ! ApplicationInitialGuess - Computes initial approximation, called by
324: ! the higher level routine FormInitialGuess().
325: !
326: ! Input Parameter:
327: ! x - local vector data
328: !
329: ! Output Parameters:
330: ! f - local vector data, f(x)
331: ! ierr - error code
332: !
333: ! Notes:
334: ! This routine uses standard Fortran-style computations over a 2-dim array.
335: !
336: subroutine ApplicationInitialGuess(x,ierr)
337: use petscksp
338: implicit none
340: ! Common blocks:
341: PetscReal lambda
342: PetscInt mx,my
343: PetscBool fd_coloring
344: common /params/ lambda,mx,my,fd_coloring
346: ! Input/output variables:
347: PetscScalar x(mx,my)
348: PetscErrorCode ierr
350: ! Local variables:
351: PetscInt i,j
352: PetscReal temp1,temp,hx,hy,one
354: ! Set parameters
356: ierr = 0
357: one = 1.0
358: hx = one/(mx-1)
359: hy = one/(my-1)
360: temp1 = lambda/(lambda + one)
362: do 20 j=1,my
363: temp = min(j-1,my-j)*hy
364: do 10 i=1,mx
365: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
366: x(i,j) = 0.0
367: else
368: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
369: endif
370: 10 continue
371: 20 continue
373: end
375: ! ---------------------------------------------------------------------
376: !
377: ! FormFunction - Evaluates nonlinear function, F(x).
378: !
379: ! Input Parameters:
380: ! snes - the SNES context
381: ! X - input vector
382: ! dummy - optional user-defined context, as set by SNESSetFunction()
383: ! (not used here)
384: !
385: ! Output Parameter:
386: ! F - vector with newly computed function
387: !
388: ! Notes:
389: ! This routine serves as a wrapper for the lower-level routine
390: ! "ApplicationFunction", where the actual computations are
391: ! done using the standard Fortran style of treating the local
392: ! vector data as a multidimensional array over the local mesh.
393: ! This routine merely accesses the local vector data via
394: ! VecGetArrayF90() and VecRestoreArrayF90().
395: !
396: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
397: use petscsnes
398: implicit none
400: ! Input/output variables:
401: SNES snes
402: Vec X,F
403: PetscErrorCode ierr
404: MatFDColoring fdcoloring
406: ! Common blocks:
407: PetscReal lambda
408: PetscInt mx,my
409: PetscBool fd_coloring
410: common /params/ lambda,mx,my,fd_coloring
412: ! Declarations for use with local arrays:
413: PetscScalar,pointer :: lx_v(:), lf_v(:)
414: PetscInt, pointer :: indices(:)
416: ! Get pointers to vector data.
417: ! - VecGetArrayF90() returns a pointer to the data array.
418: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
419: ! the array.
421: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
422: PetscCallA(VecGetArrayF90(F,lf_v,ierr))
424: ! Compute function
426: PetscCallA(ApplicationFunction(lx_v,lf_v,ierr))
428: ! Restore vectors
430: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
431: PetscCallA(VecRestoreArrayF90(F,lf_v,ierr))
433: PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
434: !
435: ! fdcoloring is in the common block and used here ONLY to test the
436: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
437: !
438: if (fd_coloring) then
439: PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
440: print*,'Indices from GetPerturbedColumnsF90'
441: write(*,1000) indices
442: 1000 format(50i4)
443: PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
444: endif
445: end
447: ! ---------------------------------------------------------------------
448: !
449: ! ApplicationFunction - Computes nonlinear function, called by
450: ! the higher level routine FormFunction().
451: !
452: ! Input Parameter:
453: ! x - local vector data
454: !
455: ! Output Parameters:
456: ! f - local vector data, f(x)
457: ! ierr - error code
458: !
459: ! Notes:
460: ! This routine uses standard Fortran-style computations over a 2-dim array.
461: !
462: subroutine ApplicationFunction(x,f,ierr)
463: use petscsnes
464: implicit none
466: ! Common blocks:
467: PetscReal lambda
468: PetscInt mx,my
469: PetscBool fd_coloring
470: common /params/ lambda,mx,my,fd_coloring
472: ! Input/output variables:
473: PetscScalar x(mx,my),f(mx,my)
474: PetscErrorCode ierr
476: ! Local variables:
477: PetscScalar two,one,hx,hy
478: PetscScalar hxdhy,hydhx,sc
479: PetscScalar u,uxx,uyy
480: PetscInt i,j
482: ierr = 0
483: one = 1.0
484: two = 2.0
485: hx = one/(mx-1)
486: hy = one/(my-1)
487: sc = hx*hy*lambda
488: hxdhy = hx/hy
489: hydhx = hy/hx
491: ! Compute function
493: do 20 j=1,my
494: do 10 i=1,mx
495: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
496: f(i,j) = x(i,j)
497: else
498: u = x(i,j)
499: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
500: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
501: f(i,j) = uxx + uyy - sc*exp(u)
502: endif
503: 10 continue
504: 20 continue
506: end
508: ! ---------------------------------------------------------------------
509: !
510: ! FormJacobian - Evaluates Jacobian matrix.
511: !
512: ! Input Parameters:
513: ! snes - the SNES context
514: ! x - input vector
515: ! dummy - optional user-defined context, as set by SNESSetJacobian()
516: ! (not used here)
517: !
518: ! Output Parameters:
519: ! jac - Jacobian matrix
520: ! jac_prec - optionally different preconditioning matrix (not used here)
521: !
522: ! Notes:
523: ! This routine serves as a wrapper for the lower-level routine
524: ! "ApplicationJacobian", where the actual computations are
525: ! done using the standard Fortran style of treating the local
526: ! vector data as a multidimensional array over the local mesh.
527: ! This routine merely accesses the local vector data via
528: ! VecGetArrayF90() and VecRestoreArrayF90().
529: !
530: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
531: use petscsnes
532: implicit none
534: ! Input/output variables:
535: SNES snes
536: Vec X
537: Mat jac,jac_prec
538: PetscErrorCode ierr
539: integer dummy
541: ! Common blocks:
542: PetscReal lambda
543: PetscInt mx,my
544: PetscBool fd_coloring
545: common /params/ lambda,mx,my,fd_coloring
547: ! Declarations for use with local array:
548: PetscScalar,pointer :: lx_v(:)
550: ! Get a pointer to vector data
552: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
554: ! Compute Jacobian entries
556: PetscCallA(ApplicationJacobian(lx_v,jac,jac_prec,ierr))
558: ! Restore vector
560: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
562: ! Assemble matrix
564: PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
565: PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
567: end
569: ! ---------------------------------------------------------------------
570: !
571: ! ApplicationJacobian - Computes Jacobian matrix, called by
572: ! the higher level routine FormJacobian().
573: !
574: ! Input Parameters:
575: ! x - local vector data
576: !
577: ! Output Parameters:
578: ! jac - Jacobian matrix
579: ! jac_prec - optionally different preconditioning matrix (not used here)
580: ! ierr - error code
581: !
582: ! Notes:
583: ! This routine uses standard Fortran-style computations over a 2-dim array.
584: !
585: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
586: use petscsnes
587: implicit none
589: ! Common blocks:
590: PetscReal lambda
591: PetscInt mx,my
592: PetscBool fd_coloring
593: common /params/ lambda,mx,my,fd_coloring
595: ! Input/output variables:
596: PetscScalar x(mx,my)
597: Mat jac,jac_prec
598: PetscErrorCode ierr
600: ! Local variables:
601: PetscInt i,j,row(1),col(5),i1,i5
602: PetscScalar two,one, hx,hy
603: PetscScalar hxdhy,hydhx,sc,v(5)
605: ! Set parameters
607: i1 = 1
608: i5 = 5
609: one = 1.0
610: two = 2.0
611: hx = one/(mx-1)
612: hy = one/(my-1)
613: sc = hx*hy
614: hxdhy = hx/hy
615: hydhx = hy/hx
617: ! Compute entries of the Jacobian matrix
618: ! - Here, we set all entries for a particular row at once.
619: ! - Note that MatSetValues() uses 0-based row and column numbers
620: ! in Fortran as well as in C.
622: do 20 j=1,my
623: row(1) = (j-1)*mx - 1
624: do 10 i=1,mx
625: row(1) = row(1) + 1
626: ! boundary points
627: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
628: PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,[one],INSERT_VALUES,ierr))
629: ! interior grid points
630: else
631: v(1) = -hxdhy
632: v(2) = -hydhx
633: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
634: v(4) = -hydhx
635: v(5) = -hxdhy
636: col(1) = row(1) - mx
637: col(2) = row(1) - 1
638: col(3) = row(1)
639: col(4) = row(1) + 1
640: col(5) = row(1) + mx
641: PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
642: endif
643: 10 continue
644: 20 continue
646: end
648: !
649: !/*TEST
650: !
651: ! build:
652: ! requires: !single
653: !
654: ! test:
655: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
656: !
657: ! test:
658: ! suffix: 2
659: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
660: !
661: ! test:
662: ! suffix: 3
663: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
664: ! filter: sort -b
665: ! filter_output: sort -b
666: !
667: ! test:
668: ! suffix: 4
669: ! args: -pc -par 6.807 -nox
670: !
671: !TEST*/