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
42: PetscCallA(VecDuplicate(x,tmp,ierr))
43: mone = -1.0
44: PetscCallA(VecWAXPY(tmp,mone,x,w,ierr))
45: PetscCallA(VecNorm(tmp,NORM_2,norm,ierr))
46: PetscCallA(VecDestroy(tmp,ierr))
47: print*, 'Norm of search step ',norm
48: changed_y = PETSC_FALSE
49: changed_w = PETSC_FALSE
50: return
51: end
53: program main
54: #include <petsc/finclude/petscdraw.h>
55: use petscsnes
56: implicit none
57: interface SNESSetJacobian
58: subroutine SNESSetJacobian1(a,b,c,d,e,z)
59: use petscsnes
60: SNES a
61: Mat b
62: Mat c
63: external d
64: MatFDColoring e
65: PetscErrorCode z
66: end subroutine
67: subroutine SNESSetJacobian2(a,b,c,d,e,z)
68: use petscsnes
69: SNES a
70: Mat b
71: Mat c
72: external d
73: integer e
74: PetscErrorCode z
75: end subroutine
76: end interface
77: !
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: ! Variable declarations
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: !
82: ! Variables:
83: ! snes - nonlinear solver
84: ! x, r - solution, residual vectors
85: ! J - Jacobian matrix
86: ! its - iterations for convergence
87: ! matrix_free - flag - 1 indicates matrix-free version
88: ! lambda - nonlinearity parameter
89: ! draw - drawing context
90: !
91: SNES snes
92: MatColoring mc
93: Vec x,r
94: PetscDraw draw
95: Mat J
96: PetscBool matrix_free,flg,fd_coloring
97: PetscErrorCode ierr
98: PetscInt its,N, mx,my,i5
99: PetscMPIInt size,rank
100: PetscReal lambda_max,lambda_min,lambda
101: MatFDColoring fdcoloring
102: ISColoring iscoloring
103: PetscBool pc
104: external postcheck
106: PetscScalar lx_v(0:1)
107: PetscOffset lx_i
109: ! Store parameters in common block
111: common /params/ lambda,mx,my,fd_coloring
113: ! Note: Any user-defined Fortran routines (such as FormJacobian)
114: ! MUST be declared as external.
116: external FormFunction,FormInitialGuess,FormJacobian
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: ! Initialize program
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: PetscCallA(PetscInitialize(ierr))
123: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
124: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
126: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
128: ! Initialize problem parameters
129: i5 = 5
130: lambda_max = 6.81
131: lambda_min = 0.0
132: lambda = 6.0
133: mx = 4
134: my = 4
135: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
136: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
137: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
138: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
139: N = mx*my
140: pc = PETSC_FALSE
141: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Create nonlinear solver context
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
149: if (pc .eqv. PETSC_TRUE) then
150: PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
151: PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
152: endif
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: ! Create vector data structures; set function evaluation routine
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
159: PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
160: PetscCallA(VecSetFromOptions(x,ierr))
161: PetscCallA(VecDuplicate(x,r,ierr))
163: ! Set function evaluation routine and vector. Whenever the nonlinear
164: ! solver needs to evaluate the nonlinear function, it will call this
165: ! routine.
166: ! - Note that the final routine argument is the user-defined
167: ! context that provides application-specific data for the
168: ! function evaluation routine.
170: PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: ! Create matrix data structure; set Jacobian evaluation routine
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: ! Create matrix. Here we only approximately preallocate storage space
177: ! for the Jacobian. See the users manual for a discussion of better
178: ! techniques for preallocating matrix memory.
180: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
181: if (.not. matrix_free) then
182: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr))
183: endif
185: !
186: ! This option will cause the Jacobian to be computed via finite differences
187: ! efficiently using a coloring of the columns of the matrix.
188: !
189: fd_coloring = .false.
190: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
191: if (fd_coloring) then
193: !
194: ! This initializes the nonzero structure of the Jacobian. This is artificial
195: ! because clearly if we had a routine to compute the Jacobian we won't need
196: ! to use finite differences.
197: !
198: PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
199: !
200: ! Color the matrix, i.e. determine groups of columns that share no common
201: ! rows. These columns in the Jacobian can all be computed simultaneously.
202: !
203: PetscCallA(MatColoringCreate(J,mc,ierr))
204: PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
205: PetscCallA(MatColoringSetFromOptions(mc,ierr))
206: PetscCallA(MatColoringApply(mc,iscoloring,ierr))
207: PetscCallA(MatColoringDestroy(mc,ierr))
208: !
209: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
210: ! to compute the actual Jacobians via finite differences.
211: !
212: PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
213: PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
214: PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
215: PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
216: !
217: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
218: ! to compute Jacobians.
219: !
220: PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
221: PetscCallA(ISColoringDestroy(iscoloring,ierr))
223: else if (.not. matrix_free) then
225: ! Set Jacobian matrix data structure and default Jacobian evaluation
226: ! routine. Whenever the nonlinear solver needs to compute the
227: ! Jacobian matrix, it will call this routine.
228: ! - Note that the final routine argument is the user-defined
229: ! context that provides application-specific data for the
230: ! Jacobian evaluation routine.
231: ! - The user can override with:
232: ! -snes_fd : default finite differencing approximation of Jacobian
233: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
234: ! (unless user explicitly sets preconditioner)
235: ! -snes_mf_operator : form preconditioning matrix as set by the user,
236: ! but use matrix-free approx for Jacobian-vector
237: ! products within Newton-Krylov method
238: !
239: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
240: endif
242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: ! Customize nonlinear solver; set runtime options
244: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
246: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
248: PetscCallA(SNESSetFromOptions(snes,ierr))
250: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: ! Evaluate initial guess; then solve nonlinear system.
252: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: ! Note: The user should initialize the vector, x, with the initial guess
255: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
256: ! to employ an initial guess of zero, the user should explicitly set
257: ! this vector to zero by calling VecSet().
259: PetscCallA(FormInitialGuess(x,ierr))
260: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
261: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
262: if (rank .eq. 0) then
263: write(6,100) its
264: endif
265: 100 format('Number of SNES iterations = ',i1)
267: ! PetscDraw contour plot of solution
269: PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
270: PetscCallA(PetscDrawSetFromOptions(draw,ierr))
272: PetscCallA(VecGetArrayRead(x,lx_v,lx_i,ierr))
273: PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr))
274: PetscCallA(VecRestoreArrayRead(x,lx_v,lx_i,ierr))
276: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: ! Free work space. All PETSc objects should be destroyed when they
278: ! are no longer needed.
279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
282: if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))
284: PetscCallA(VecDestroy(x,ierr))
285: PetscCallA(VecDestroy(r,ierr))
286: PetscCallA(SNESDestroy(snes,ierr))
287: PetscCallA(PetscDrawDestroy(draw,ierr))
288: PetscCallA(PetscFinalize(ierr))
289: end
291: ! ---------------------------------------------------------------------
292: !
293: ! FormInitialGuess - Forms initial approximation.
294: !
295: ! Input Parameter:
296: ! X - vector
297: !
298: ! Output Parameters:
299: ! X - vector
300: ! ierr - error code
301: !
302: ! Notes:
303: ! This routine serves as a wrapper for the lower-level routine
304: ! "ApplicationInitialGuess", where the actual computations are
305: ! done using the standard Fortran style of treating the local
306: ! vector data as a multidimensional array over the local mesh.
307: ! This routine merely accesses the local vector data via
308: ! VecGetArray() and VecRestoreArray().
309: !
310: subroutine FormInitialGuess(X,ierr)
311: use petscsnes
312: implicit none
314: ! Input/output variables:
315: Vec X
316: PetscErrorCode ierr
318: ! Declarations for use with local arrays:
319: PetscScalar lx_v(0:1)
320: PetscOffset lx_i
322: 0
324: ! Get a pointer to vector data.
325: ! - For default PETSc vectors, VecGetArray() returns a pointer to
326: ! the data array. Otherwise, the routine is implementation dependent.
327: ! - You MUST call VecRestoreArray() when you no longer need access to
328: ! the array.
329: ! - Note that the Fortran interface to VecGetArray() differs from the
330: ! C version. See the users manual for details.
332: PetscCallA(VecGetArray(X,lx_v,lx_i,ierr))
334: ! Compute initial guess
336: PetscCallA(ApplicationInitialGuess(lx_v(lx_i),ierr))
338: ! Restore vector
340: PetscCallA(VecRestoreArray(X,lx_v,lx_i,ierr))
342: return
343: end
345: ! ApplicationInitialGuess - Computes initial approximation, called by
346: ! the higher level routine FormInitialGuess().
347: !
348: ! Input Parameter:
349: ! x - local vector data
350: !
351: ! Output Parameters:
352: ! f - local vector data, f(x)
353: ! ierr - error code
354: !
355: ! Notes:
356: ! This routine uses standard Fortran-style computations over a 2-dim array.
357: !
358: subroutine ApplicationInitialGuess(x,ierr)
359: use petscksp
360: implicit none
362: ! Common blocks:
363: PetscReal lambda
364: PetscInt mx,my
365: PetscBool fd_coloring
366: common /params/ lambda,mx,my,fd_coloring
368: ! Input/output variables:
369: PetscScalar x(mx,my)
370: PetscErrorCode ierr
372: ! Local variables:
373: PetscInt i,j
374: PetscReal temp1,temp,hx,hy,one
376: ! Set parameters
378: 0
379: one = 1.0
380: hx = one/(mx-1)
381: hy = one/(my-1)
382: temp1 = lambda/(lambda + one)
384: do 20 j=1,my
385: temp = min(j-1,my-j)*hy
386: do 10 i=1,mx
387: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
388: x(i,j) = 0.0
389: else
390: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
391: endif
392: 10 continue
393: 20 continue
395: return
396: end
398: ! ---------------------------------------------------------------------
399: !
400: ! FormFunction - Evaluates nonlinear function, F(x).
401: !
402: ! Input Parameters:
403: ! snes - the SNES context
404: ! X - input vector
405: ! dummy - optional user-defined context, as set by SNESSetFunction()
406: ! (not used here)
407: !
408: ! Output Parameter:
409: ! F - vector with newly computed function
410: !
411: ! Notes:
412: ! This routine serves as a wrapper for the lower-level routine
413: ! "ApplicationFunction", where the actual computations are
414: ! done using the standard Fortran style of treating the local
415: ! vector data as a multidimensional array over the local mesh.
416: ! This routine merely accesses the local vector data via
417: ! VecGetArray() and VecRestoreArray().
418: !
419: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
420: use petscsnes
421: implicit none
423: ! Input/output variables:
424: SNES snes
425: Vec X,F
426: PetscErrorCode ierr
427: MatFDColoring fdcoloring
429: ! Common blocks:
430: PetscReal lambda
431: PetscInt mx,my
432: PetscBool fd_coloring
433: common /params/ lambda,mx,my,fd_coloring
435: ! Declarations for use with local arrays:
436: PetscScalar lx_v(0:1),lf_v(0:1)
437: PetscOffset lx_i,lf_i
439: PetscInt, pointer :: indices(:)
441: ! Get pointers to vector data.
442: ! - For default PETSc vectors, VecGetArray() returns a pointer to
443: ! the data array. Otherwise, the routine is implementation dependent.
444: ! - You MUST call VecRestoreArray() when you no longer need access to
445: ! the array.
446: ! - Note that the Fortran interface to VecGetArray() differs from the
447: ! C version. See the Fortran chapter of the users manual for details.
449: PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr))
450: PetscCallA(VecGetArray(F,lf_v,lf_i,ierr))
452: ! Compute function
454: PetscCallA(ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr))
456: ! Restore vectors
458: PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr))
459: PetscCallA(VecRestoreArray(F,lf_v,lf_i,ierr))
461: PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
462: !
463: ! fdcoloring is in the common block and used here ONLY to test the
464: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
465: !
466: if (fd_coloring) then
467: PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
468: print*,'Indices from GetPerturbedColumnsF90'
469: write(*,1000) indices
470: 1000 format(50i4)
471: PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
472: endif
473: return
474: end
476: ! ---------------------------------------------------------------------
477: !
478: ! ApplicationFunction - Computes nonlinear function, called by
479: ! the higher level routine FormFunction().
480: !
481: ! Input Parameter:
482: ! x - local vector data
483: !
484: ! Output Parameters:
485: ! f - local vector data, f(x)
486: ! ierr - error code
487: !
488: ! Notes:
489: ! This routine uses standard Fortran-style computations over a 2-dim array.
490: !
491: subroutine ApplicationFunction(x,f,ierr)
492: use petscsnes
493: implicit none
495: ! Common blocks:
496: PetscReal lambda
497: PetscInt mx,my
498: PetscBool fd_coloring
499: common /params/ lambda,mx,my,fd_coloring
501: ! Input/output variables:
502: PetscScalar x(mx,my),f(mx,my)
503: PetscErrorCode ierr
505: ! Local variables:
506: PetscScalar two,one,hx,hy
507: PetscScalar hxdhy,hydhx,sc
508: PetscScalar u,uxx,uyy
509: PetscInt i,j
511: 0
512: one = 1.0
513: two = 2.0
514: hx = one/(mx-1)
515: hy = one/(my-1)
516: sc = hx*hy*lambda
517: hxdhy = hx/hy
518: hydhx = hy/hx
520: ! Compute function
522: do 20 j=1,my
523: do 10 i=1,mx
524: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
525: f(i,j) = x(i,j)
526: else
527: u = x(i,j)
528: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
529: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
530: f(i,j) = uxx + uyy - sc*exp(u)
531: endif
532: 10 continue
533: 20 continue
535: return
536: end
538: ! ---------------------------------------------------------------------
539: !
540: ! FormJacobian - Evaluates Jacobian matrix.
541: !
542: ! Input Parameters:
543: ! snes - the SNES context
544: ! x - input vector
545: ! dummy - optional user-defined context, as set by SNESSetJacobian()
546: ! (not used here)
547: !
548: ! Output Parameters:
549: ! jac - Jacobian matrix
550: ! jac_prec - optionally different preconditioning matrix (not used here)
551: ! flag - flag indicating matrix structure
552: !
553: ! Notes:
554: ! This routine serves as a wrapper for the lower-level routine
555: ! "ApplicationJacobian", where the actual computations are
556: ! done using the standard Fortran style of treating the local
557: ! vector data as a multidimensional array over the local mesh.
558: ! This routine merely accesses the local vector data via
559: ! VecGetArray() and VecRestoreArray().
560: !
561: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
562: use petscsnes
563: implicit none
565: ! Input/output variables:
566: SNES snes
567: Vec X
568: Mat jac,jac_prec
569: PetscErrorCode ierr
570: integer dummy
572: ! Common blocks:
573: PetscReal lambda
574: PetscInt mx,my
575: PetscBool fd_coloring
576: common /params/ lambda,mx,my,fd_coloring
578: ! Declarations for use with local array:
579: PetscScalar lx_v(0:1)
580: PetscOffset lx_i
582: ! Get a pointer to vector data
584: PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr))
586: ! Compute Jacobian entries
588: PetscCallA(ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr))
590: ! Restore vector
592: PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr))
594: ! Assemble matrix
596: PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
597: PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
599: return
600: end
602: ! ---------------------------------------------------------------------
603: !
604: ! ApplicationJacobian - Computes Jacobian matrix, called by
605: ! the higher level routine FormJacobian().
606: !
607: ! Input Parameters:
608: ! x - local vector data
609: !
610: ! Output Parameters:
611: ! jac - Jacobian matrix
612: ! jac_prec - optionally different preconditioning matrix (not used here)
613: ! ierr - error code
614: !
615: ! Notes:
616: ! This routine uses standard Fortran-style computations over a 2-dim array.
617: !
618: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
619: use petscsnes
620: implicit none
622: ! Common blocks:
623: PetscReal lambda
624: PetscInt mx,my
625: PetscBool fd_coloring
626: common /params/ lambda,mx,my,fd_coloring
628: ! Input/output variables:
629: PetscScalar x(mx,my)
630: Mat jac,jac_prec
631: PetscErrorCode ierr
633: ! Local variables:
634: PetscInt i,j,row(1),col(5),i1,i5
635: PetscScalar two,one, hx,hy
636: PetscScalar hxdhy,hydhx,sc,v(5)
638: ! Set parameters
640: i1 = 1
641: i5 = 5
642: one = 1.0
643: two = 2.0
644: hx = one/(mx-1)
645: hy = one/(my-1)
646: sc = hx*hy
647: hxdhy = hx/hy
648: hydhx = hy/hx
650: ! Compute entries of the Jacobian matrix
651: ! - Here, we set all entries for a particular row at once.
652: ! - Note that MatSetValues() uses 0-based row and column numbers
653: ! in Fortran as well as in C.
655: do 20 j=1,my
656: row(1) = (j-1)*mx - 1
657: do 10 i=1,mx
658: row(1) = row(1) + 1
659: ! boundary points
660: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
661: PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr))
662: ! interior grid points
663: else
664: v(1) = -hxdhy
665: v(2) = -hydhx
666: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
667: v(4) = -hydhx
668: v(5) = -hxdhy
669: col(1) = row(1) - mx
670: col(2) = row(1) - 1
671: col(3) = row(1)
672: col(4) = row(1) + 1
673: col(5) = row(1) + mx
674: PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
675: endif
676: 10 continue
677: 20 continue
679: return
680: end
682: !
683: !/*TEST
684: !
685: ! build:
686: ! requires: !single
687: !
688: ! test:
689: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
690: !
691: ! test:
692: ! suffix: 2
693: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
694: !
695: ! test:
696: ! suffix: 3
697: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
698: ! filter: sort -b
699: ! filter_output: sort -b
700: !
701: ! test:
702: ! suffix: 4
703: ! args: -pc -par 6.807 -nox
704: !
705: !TEST*/