Actual source code: ex1f.F90
petsc-3.13.4 2020-08-01
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: !
10: !!/*T
11: ! Concepts: SNES^sequential Bratu example
12: ! Processors: 1
13: !T*/
16: !
17: ! --------------------------------------------------------------------------
18: !
19: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
20: ! the partial differential equation
21: !
22: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
23: !
24: ! with boundary conditions
25: !
26: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
27: !
28: ! A finite difference approximation with the usual 5-point stencil
29: ! is used to discretize the boundary value problem to obtain a nonlinear
30: ! system of equations.
31: !
32: ! The parallel version of this code is snes/tutorials/ex5f.F
33: !
34: ! --------------------------------------------------------------------------
35: subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr)
36: #include <petsc/finclude/petscsnes.h>
37: use petscsnes
38: implicit none
39: SNES snes
40: PetscReal norm
41: Vec tmp,x,y,w
42: PetscBool changed_w,changed_y
43: PetscErrorCode ierr
44: PetscInt ctx
45: PetscScalar mone
47: call VecDuplicate(x,tmp,ierr)
48: mone = -1.0
49: call VecWAXPY(tmp,mone,x,w,ierr)
50: call VecNorm(tmp,NORM_2,norm,ierr)
51: call VecDestroy(tmp,ierr)
52: print*, 'Norm of search step ',norm
53: changed_y = PETSC_FALSE
54: changed_w = PETSC_FALSE
55: return
56: end
58: program main
59: #include <petsc/finclude/petscdraw.h>
60: use petscsnes
61: implicit none
62: interface SNESSetJacobian
63: subroutine SNESSetJacobian1(a,b,c,d,e,z)
64: use petscsnes
65: SNES a
66: Mat b
67: Mat c
68: external d
69: MatFDColoring e
70: PetscErrorCode z
71: end subroutine
72: subroutine SNESSetJacobian2(a,b,c,d,e,z)
73: use petscsnes
74: SNES a
75: Mat b
76: Mat c
77: external d
78: integer e
79: PetscErrorCode z
80: end subroutine
81: end interface
82: !
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Variable declarations
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: !
87: ! Variables:
88: ! snes - nonlinear solver
89: ! x, r - solution, residual vectors
90: ! J - Jacobian matrix
91: ! its - iterations for convergence
92: ! matrix_free - flag - 1 indicates matrix-free version
93: ! lambda - nonlinearity parameter
94: ! draw - drawing context
95: !
96: SNES snes
97: MatColoring mc
98: Vec x,r
99: PetscDraw draw
100: Mat J
101: PetscBool matrix_free,flg,fd_coloring
102: PetscErrorCode ierr
103: PetscInt its,N, mx,my,i5
104: PetscMPIInt size,rank
105: PetscReal lambda_max,lambda_min,lambda
106: MatFDColoring fdcoloring
107: ISColoring iscoloring
108: PetscBool pc
109: external postcheck
111: PetscScalar lx_v(0:1)
112: PetscOffset lx_i
114: ! Store parameters in common block
116: common /params/ lambda,mx,my,fd_coloring
118: ! Note: Any user-defined Fortran routines (such as FormJacobian)
119: ! MUST be declared as external.
121: external FormFunction,FormInitialGuess,FormJacobian
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! Initialize program
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
128: if (ierr .ne. 0) then
129: print*,'Unable to initialize PETSc'
130: stop
131: endif
132: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
133: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
135: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
137: ! Initialize problem parameters
138: i5 = 5
139: lambda_max = 6.81
140: lambda_min = 0.0
141: lambda = 6.0
142: mx = 4
143: my = 4
144: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
145: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
146: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
147: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
148: N = mx*my
149: pc = PETSC_FALSE
150: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr);
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! Create nonlinear solver context
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
158: if (pc .eqv. PETSC_TRUE) then
159: call SNESSetType(snes,SNESNEWTONTR,ierr)
160: call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
161: endif
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Create vector data structures; set function evaluation routine
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: call VecCreate(PETSC_COMM_WORLD,x,ierr)
168: call VecSetSizes(x,PETSC_DECIDE,N,ierr)
169: call VecSetFromOptions(x,ierr)
170: call VecDuplicate(x,r,ierr)
172: ! Set function evaluation routine and vector. Whenever the nonlinear
173: ! solver needs to evaluate the nonlinear function, it will call this
174: ! routine.
175: ! - Note that the final routine argument is the user-defined
176: ! context that provides Section 1.5 Writing Application Codes with PETSc-specific data for the
177: ! function evaluation routine.
179: call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)
181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: ! Create matrix data structure; set Jacobian evaluation routine
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! Create matrix. Here we only approximately preallocate storage space
186: ! for the Jacobian. See the users manual for a discussion of better
187: ! techniques for preallocating matrix memory.
189: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
190: if (.not. matrix_free) then
191: call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
192: endif
194: !
195: ! This option will cause the Jacobian to be computed via finite differences
196: ! efficiently using a coloring of the columns of the matrix.
197: !
198: fd_coloring = .false.
199: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
200: if (fd_coloring) then
202: !
203: ! This initializes the nonzero structure of the Jacobian. This is artificial
204: ! because clearly if we had a routine to compute the Jacobian we won't need
205: ! to use finite differences.
206: !
207: call FormJacobian(snes,x,J,J,0,ierr)
208: !
209: ! Color the matrix, i.e. determine groups of columns that share no common
210: ! rows. These columns in the Jacobian can all be computed simulataneously.
211: !
212: call MatColoringCreate(J,mc,ierr)
213: call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
214: call MatColoringSetFromOptions(mc,ierr)
215: call MatColoringApply(mc,iscoloring,ierr)
216: call MatColoringDestroy(mc,ierr)
217: !
218: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
219: ! to compute the actual Jacobians via finite differences.
220: !
221: call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
222: call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
223: call MatFDColoringSetFromOptions(fdcoloring,ierr)
224: call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
225: !
226: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
227: ! to compute Jacobians.
228: !
229: call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
230: call ISColoringDestroy(iscoloring,ierr)
232: else if (.not. matrix_free) then
234: ! Set Jacobian matrix data structure and default Jacobian evaluation
235: ! routine. Whenever the nonlinear solver needs to compute the
236: ! Jacobian matrix, it will call this routine.
237: ! - Note that the final routine argument is the user-defined
238: ! context that provides Section 1.5 Writing Application Codes with PETSc-specific data for the
239: ! Jacobian evaluation routine.
240: ! - The user can override with:
241: ! -snes_fd : default finite differencing approximation of Jacobian
242: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
243: ! (unless user explicitly sets preconditioner)
244: ! -snes_mf_operator : form preconditioning matrix as set by the user,
245: ! but use matrix-free approx for Jacobian-vector
246: ! products within Newton-Krylov method
247: !
248: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
249: endif
251: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: ! Customize nonlinear solver; set runtime options
253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
257: call SNESSetFromOptions(snes,ierr)
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: ! Evaluate initial guess; then solve nonlinear system.
261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263: ! Note: The user should initialize the vector, x, with the initial guess
264: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
265: ! to employ an initial guess of zero, the user should explicitly set
266: ! this vector to zero by calling VecSet().
268: call FormInitialGuess(x,ierr)
269: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
270: call SNESGetIterationNumber(snes,its,ierr);
271: if (rank .eq. 0) then
272: write(6,100) its
273: endif
274: 100 format('Number of SNES iterations = ',i1)
276: ! PetscDraw contour plot of solution
278: call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
279: call PetscDrawSetFromOptions(draw,ierr)
281: call VecGetArrayRead(x,lx_v,lx_i,ierr)
282: call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
283: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
285: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: ! Free work space. All PETSc objects should be destroyed when they
287: ! are no longer needed.
288: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
290: if (.not. matrix_free) call MatDestroy(J,ierr)
291: if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)
293: call VecDestroy(x,ierr)
294: call VecDestroy(r,ierr)
295: call SNESDestroy(snes,ierr)
296: call PetscDrawDestroy(draw,ierr)
297: call PetscFinalize(ierr)
298: end
300: ! ---------------------------------------------------------------------
301: !
302: ! FormInitialGuess - Forms initial approximation.
303: !
304: ! Input Parameter:
305: ! X - vector
306: !
307: ! Output Parameters:
308: ! X - vector
309: ! ierr - error code
310: !
311: ! Notes:
312: ! This routine serves as a wrapper for the lower-level routine
313: ! "ApplicationInitialGuess", where the actual computations are
314: ! done using the standard Fortran style of treating the local
315: ! vector data as a multidimensional array over the local mesh.
316: ! This routine merely accesses the local vector data via
317: ! VecGetArray() and VecRestoreArray().
318: !
319: subroutine FormInitialGuess(X,ierr)
320: use petscsnes
321: implicit none
323: ! Input/output variables:
324: Vec X
325: PetscErrorCode ierr
327: ! Declarations for use with local arrays:
328: PetscScalar lx_v(0:1)
329: PetscOffset lx_i
331: 0
333: ! Get a pointer to vector data.
334: ! - For default PETSc vectors, VecGetArray() returns a pointer to
335: ! the data array. Otherwise, the routine is implementation dependent.
336: ! - You MUST call VecRestoreArray() when you no longer need access to
337: ! the array.
338: ! - Note that the Fortran interface to VecGetArray() differs from the
339: ! C version. See the users manual for details.
341: call VecGetArray(X,lx_v,lx_i,ierr)
343: ! Compute initial guess
345: call ApplicationInitialGuess(lx_v(lx_i),ierr)
347: ! Restore vector
349: call VecRestoreArray(X,lx_v,lx_i,ierr)
351: return
352: end
354: ! ---------------------------------------------------------------------
355: !
356: ! ApplicationInitialGuess - Computes initial approximation, called by
357: ! the higher level routine FormInitialGuess().
358: !
359: ! Input Parameter:
360: ! x - local vector data
361: !
362: ! Output Parameters:
363: ! f - local vector data, f(x)
364: ! ierr - error code
365: !
366: ! Notes:
367: ! This routine uses standard Fortran-style computations over a 2-dim array.
368: !
369: subroutine ApplicationInitialGuess(x,ierr)
370: use petscksp
371: implicit none
373: ! Common blocks:
374: PetscReal lambda
375: PetscInt mx,my
376: PetscBool fd_coloring
377: common /params/ lambda,mx,my,fd_coloring
379: ! Input/output variables:
380: PetscScalar x(mx,my)
381: PetscErrorCode ierr
383: ! Local variables:
384: PetscInt i,j
385: PetscReal temp1,temp,hx,hy,one
387: ! Set parameters
389: 0
390: one = 1.0
391: hx = one/(mx-1)
392: hy = one/(my-1)
393: temp1 = lambda/(lambda + one)
395: do 20 j=1,my
396: temp = min(j-1,my-j)*hy
397: do 10 i=1,mx
398: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
399: x(i,j) = 0.0
400: else
401: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
402: endif
403: 10 continue
404: 20 continue
406: return
407: end
409: ! ---------------------------------------------------------------------
410: !
411: ! FormFunction - Evaluates nonlinear function, F(x).
412: !
413: ! Input Parameters:
414: ! snes - the SNES context
415: ! X - input vector
416: ! dummy - optional user-defined context, as set by SNESSetFunction()
417: ! (not used here)
418: !
419: ! Output Parameter:
420: ! F - vector with newly computed function
421: !
422: ! Notes:
423: ! This routine serves as a wrapper for the lower-level routine
424: ! "ApplicationFunction", where the actual computations are
425: ! done using the standard Fortran style of treating the local
426: ! vector data as a multidimensional array over the local mesh.
427: ! This routine merely accesses the local vector data via
428: ! VecGetArray() and VecRestoreArray().
429: !
430: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
431: use petscsnes
432: implicit none
434: ! Input/output variables:
435: SNES snes
436: Vec X,F
437: PetscErrorCode ierr
438: MatFDColoring fdcoloring
440: ! Common blocks:
441: PetscReal lambda
442: PetscInt mx,my
443: PetscBool fd_coloring
444: common /params/ lambda,mx,my,fd_coloring
446: ! Declarations for use with local arrays:
447: PetscScalar lx_v(0:1),lf_v(0:1)
448: PetscOffset lx_i,lf_i
450: PetscInt, pointer :: indices(:)
452: ! Get pointers to vector data.
453: ! - For default PETSc vectors, VecGetArray() returns a pointer to
454: ! the data array. Otherwise, the routine is implementation dependent.
455: ! - You MUST call VecRestoreArray() when you no longer need access to
456: ! the array.
457: ! - Note that the Fortran interface to VecGetArray() differs from the
458: ! C version. See the Fortran chapter of the users manual for details.
460: call VecGetArrayRead(X,lx_v,lx_i,ierr)
461: call VecGetArray(F,lf_v,lf_i,ierr)
463: ! Compute function
465: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)
467: ! Restore vectors
469: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
470: call VecRestoreArray(F,lf_v,lf_i,ierr)
472: call PetscLogFlops(11.0d0*mx*my,ierr)
473: !
474: ! fdcoloring is in the common block and used here ONLY to test the
475: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
476: !
477: if (fd_coloring) then
478: call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
479: print*,'Indices from GetPerturbedColumnsF90'
480: write(*,1000) indices
481: 1000 format(50i4)
482: call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
483: endif
484: return
485: end
487: ! ---------------------------------------------------------------------
488: !
489: ! ApplicationFunction - Computes nonlinear function, called by
490: ! the higher level routine FormFunction().
491: !
492: ! Input Parameter:
493: ! x - local vector data
494: !
495: ! Output Parameters:
496: ! f - local vector data, f(x)
497: ! ierr - error code
498: !
499: ! Notes:
500: ! This routine uses standard Fortran-style computations over a 2-dim array.
501: !
502: subroutine ApplicationFunction(x,f,ierr)
503: use petscsnes
504: implicit none
506: ! Common blocks:
507: PetscReal lambda
508: PetscInt mx,my
509: PetscBool fd_coloring
510: common /params/ lambda,mx,my,fd_coloring
512: ! Input/output variables:
513: PetscScalar x(mx,my),f(mx,my)
514: PetscErrorCode ierr
516: ! Local variables:
517: PetscScalar two,one,hx,hy
518: PetscScalar hxdhy,hydhx,sc
519: PetscScalar u,uxx,uyy
520: PetscInt i,j
522: 0
523: one = 1.0
524: two = 2.0
525: hx = one/(mx-1)
526: hy = one/(my-1)
527: sc = hx*hy*lambda
528: hxdhy = hx/hy
529: hydhx = hy/hx
531: ! Compute function
533: do 20 j=1,my
534: do 10 i=1,mx
535: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
536: f(i,j) = x(i,j)
537: else
538: u = x(i,j)
539: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
540: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
541: f(i,j) = uxx + uyy - sc*exp(u)
542: endif
543: 10 continue
544: 20 continue
546: return
547: end
549: ! ---------------------------------------------------------------------
550: !
551: ! FormJacobian - Evaluates Jacobian matrix.
552: !
553: ! Input Parameters:
554: ! snes - the SNES context
555: ! x - input vector
556: ! dummy - optional user-defined context, as set by SNESSetJacobian()
557: ! (not used here)
558: !
559: ! Output Parameters:
560: ! jac - Jacobian matrix
561: ! jac_prec - optionally different preconditioning matrix (not used here)
562: ! flag - flag indicating matrix structure
563: !
564: ! Notes:
565: ! This routine serves as a wrapper for the lower-level routine
566: ! "ApplicationJacobian", where the actual computations are
567: ! done using the standard Fortran style of treating the local
568: ! vector data as a multidimensional array over the local mesh.
569: ! This routine merely accesses the local vector data via
570: ! VecGetArray() and VecRestoreArray().
571: !
572: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
573: use petscsnes
574: implicit none
576: ! Input/output variables:
577: SNES snes
578: Vec X
579: Mat jac,jac_prec
580: PetscErrorCode ierr
581: integer dummy
583: ! Common blocks:
584: PetscReal lambda
585: PetscInt mx,my
586: PetscBool fd_coloring
587: common /params/ lambda,mx,my,fd_coloring
589: ! Declarations for use with local array:
590: PetscScalar lx_v(0:1)
591: PetscOffset lx_i
593: ! Get a pointer to vector data
595: call VecGetArrayRead(X,lx_v,lx_i,ierr)
597: ! Compute Jacobian entries
599: call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)
601: ! Restore vector
603: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
605: ! Assemble matrix
607: call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
608: call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
611: return
612: end
614: ! ---------------------------------------------------------------------
615: !
616: ! ApplicationJacobian - Computes Jacobian matrix, called by
617: ! the higher level routine FormJacobian().
618: !
619: ! Input Parameters:
620: ! x - local vector data
621: !
622: ! Output Parameters:
623: ! jac - Jacobian matrix
624: ! jac_prec - optionally different preconditioning matrix (not used here)
625: ! ierr - error code
626: !
627: ! Notes:
628: ! This routine uses standard Fortran-style computations over a 2-dim array.
629: !
630: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
631: use petscsnes
632: implicit none
634: ! Common blocks:
635: PetscReal lambda
636: PetscInt mx,my
637: PetscBool fd_coloring
638: common /params/ lambda,mx,my,fd_coloring
640: ! Input/output variables:
641: PetscScalar x(mx,my)
642: Mat jac,jac_prec
643: PetscErrorCode ierr
645: ! Local variables:
646: PetscInt i,j,row(1),col(5),i1,i5
647: PetscScalar two,one, hx,hy
648: PetscScalar hxdhy,hydhx,sc,v(5)
650: ! Set parameters
652: i1 = 1
653: i5 = 5
654: one = 1.0
655: two = 2.0
656: hx = one/(mx-1)
657: hy = one/(my-1)
658: sc = hx*hy
659: hxdhy = hx/hy
660: hydhx = hy/hx
662: ! Compute entries of the Jacobian matrix
663: ! - Here, we set all entries for a particular row at once.
664: ! - Note that MatSetValues() uses 0-based row and column numbers
665: ! in Fortran as well as in C.
667: do 20 j=1,my
668: row(1) = (j-1)*mx - 1
669: do 10 i=1,mx
670: row(1) = row(1) + 1
671: ! boundary points
672: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
673: call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
674: ! interior grid points
675: else
676: v(1) = -hxdhy
677: v(2) = -hydhx
678: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
679: v(4) = -hydhx
680: v(5) = -hxdhy
681: col(1) = row(1) - mx
682: col(2) = row(1) - 1
683: col(3) = row(1)
684: col(4) = row(1) + 1
685: col(5) = row(1) + mx
686: call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
687: endif
688: 10 continue
689: 20 continue
691: return
692: end
694: !
695: !/*TEST
696: !
697: ! build:
698: ! requires: !single
699: !
700: ! test:
701: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
702: !
703: ! test:
704: ! suffix: 2
705: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
706: !
707: ! test:
708: ! suffix: 3
709: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
710: ! filter: sort -b
711: ! filter_output: sort -b
712: !
713: ! test:
714: ! suffix: 4
715: ! args: -pc -par 6.807 -nox
716: !
717: !TEST*/