Actual source code: ex1f.F90
petsc-3.10.3 2018-12-18
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/examples/tutorials/ex5f.F
33: !
34: ! --------------------------------------------------------------------------
36: program main
37: #include <petsc/finclude/petscdraw.h>
38: #include <petsc/finclude/petscsnes.h>
39: use petscsnes
40: implicit none
41: #if defined(PETSC_USING_F90) && !defined(PETSC_USE_FORTRANKIND)
42: external SNESCOMPUTEJACOBIANDEFAULTCOLOR
43: #endif
44: !
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: ! Variable declarations
47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: !
49: ! Variables:
50: ! snes - nonlinear solver
51: ! x, r - solution, residual vectors
52: ! J - Jacobian matrix
53: ! its - iterations for convergence
54: ! matrix_free - flag - 1 indicates matrix-free version
55: ! lambda - nonlinearity parameter
56: ! draw - drawing context
57: !
58: SNES snes
59: MatColoring mc
60: Vec x,r
61: PetscDraw draw
62: Mat J
63: PetscBool matrix_free,flg,fd_coloring
64: PetscErrorCode ierr
65: PetscInt its,N, mx,my,i5
66: PetscMPIInt size,rank
67: PetscReal lambda_max,lambda_min,lambda
68: MatFDColoring fdcoloring
69: ISColoring iscoloring
71: PetscScalar lx_v(0:1)
72: PetscOffset lx_i
74: ! Store parameters in common block
76: common /params/ lambda,mx,my,fd_coloring
78: ! Note: Any user-defined Fortran routines (such as FormJacobian)
79: ! MUST be declared as external.
81: external FormFunction,FormInitialGuess,FormJacobian
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Initialize program
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
88: if (ierr .ne. 0) then
89: print*,'Unable to initialize PETSc'
90: stop
91: endif
92: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
93: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
95: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only'); endif
97: ! Initialize problem parameters
98: i5 = 5
99: lambda_max = 6.81
100: lambda_min = 0.0
101: lambda = 6.0
102: mx = 4
103: my = 4
104: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
105: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
106: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
107: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,1,'Lambda out of range '); endif
108: N = mx*my
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Create nonlinear solver context
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: ! Create vector data structures; set function evaluation routine
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: call VecCreate(PETSC_COMM_WORLD,x,ierr)
121: call VecSetSizes(x,PETSC_DECIDE,N,ierr)
122: call VecSetFromOptions(x,ierr)
123: call VecDuplicate(x,r,ierr)
125: ! Set function evaluation routine and vector. Whenever the nonlinear
126: ! solver needs to evaluate the nonlinear function, it will call this
127: ! routine.
128: ! - Note that the final routine argument is the user-defined
129: ! context that provides application-specific data for the
130: ! function evaluation routine.
132: call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: ! Create matrix data structure; set Jacobian evaluation routine
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: ! Create matrix. Here we only approximately preallocate storage space
139: ! for the Jacobian. See the users manual for a discussion of better
140: ! techniques for preallocating matrix memory.
142: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
143: if (.not. matrix_free) then
144: call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
145: endif
147: !
148: ! This option will cause the Jacobian to be computed via finite differences
149: ! efficiently using a coloring of the columns of the matrix.
150: !
151: fd_coloring = .false.
152: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
153: if (fd_coloring) then
155: !
156: ! This initializes the nonzero structure of the Jacobian. This is artificial
157: ! because clearly if we had a routine to compute the Jacobian we won't need
158: ! to use finite differences.
159: !
160: call FormJacobian(snes,x,J,J,0,ierr)
161: !
162: ! Color the matrix, i.e. determine groups of columns that share no common
163: ! rows. These columns in the Jacobian can all be computed simulataneously.
164: !
165: call MatColoringCreate(J,mc,ierr)
166: call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
167: call MatColoringSetFromOptions(mc,ierr)
168: call MatColoringApply(mc,iscoloring,ierr)
169: call MatColoringDestroy(mc,ierr)
170: !
171: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
172: ! to compute the actual Jacobians via finite differences.
173: !
174: call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
175: call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
176: call MatFDColoringSetFromOptions(fdcoloring,ierr)
177: call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
178: !
179: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
180: ! to compute Jacobians.
181: !
182: call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
183: call ISColoringDestroy(iscoloring,ierr)
185: else if (.not. matrix_free) then
187: ! Set Jacobian matrix data structure and default Jacobian evaluation
188: ! routine. Whenever the nonlinear solver needs to compute the
189: ! Jacobian matrix, it will call this routine.
190: ! - Note that the final routine argument is the user-defined
191: ! context that provides application-specific data for the
192: ! Jacobian evaluation routine.
193: ! - The user can override with:
194: ! -snes_fd : default finite differencing approximation of Jacobian
195: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
196: ! (unless user explicitly sets preconditioner)
197: ! -snes_mf_operator : form preconditioning matrix as set by the user,
198: ! but use matrix-free approx for Jacobian-vector
199: ! products within Newton-Krylov method
200: !
201: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
202: endif
204: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: ! Customize nonlinear solver; set runtime options
206: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
210: call SNESSetFromOptions(snes,ierr)
212: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: ! Evaluate initial guess; then solve nonlinear system.
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: ! Note: The user should initialize the vector, x, with the initial guess
217: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
218: ! to employ an initial guess of zero, the user should explicitly set
219: ! this vector to zero by calling VecSet().
221: call FormInitialGuess(x,ierr)
222: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
223: call SNESGetIterationNumber(snes,its,ierr);
224: if (rank .eq. 0) then
225: write(6,100) its
226: endif
227: 100 format('Number of SNES iterations = ',i1)
229: ! PetscDraw contour plot of solution
231: call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
232: call PetscDrawSetFromOptions(draw,ierr)
234: call VecGetArrayRead(x,lx_v,lx_i,ierr)
235: call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
236: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
238: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: ! Free work space. All PETSc objects should be destroyed when they
240: ! are no longer needed.
241: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: if (.not. matrix_free) call MatDestroy(J,ierr)
244: if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)
246: call VecDestroy(x,ierr)
247: call VecDestroy(r,ierr)
248: call SNESDestroy(snes,ierr)
249: call PetscDrawDestroy(draw,ierr)
250: call PetscFinalize(ierr)
251: end
253: ! ---------------------------------------------------------------------
254: !
255: ! FormInitialGuess - Forms initial approximation.
256: !
257: ! Input Parameter:
258: ! X - vector
259: !
260: ! Output Parameters:
261: ! X - vector
262: ! ierr - error code
263: !
264: ! Notes:
265: ! This routine serves as a wrapper for the lower-level routine
266: ! "ApplicationInitialGuess", where the actual computations are
267: ! done using the standard Fortran style of treating the local
268: ! vector data as a multidimensional array over the local mesh.
269: ! This routine merely accesses the local vector data via
270: ! VecGetArray() and VecRestoreArray().
271: !
272: subroutine FormInitialGuess(X,ierr)
273: use petscsnes
274: implicit none
276: ! Input/output variables:
277: Vec X
278: PetscErrorCode ierr
280: ! Declarations for use with local arrays:
281: PetscScalar lx_v(0:1)
282: PetscOffset lx_i
284: 0
286: ! Get a pointer to vector data.
287: ! - For default PETSc vectors, VecGetArray() returns a pointer to
288: ! the data array. Otherwise, the routine is implementation dependent.
289: ! - You MUST call VecRestoreArray() when you no longer need access to
290: ! the array.
291: ! - Note that the Fortran interface to VecGetArray() differs from the
292: ! C version. See the users manual for details.
294: call VecGetArray(X,lx_v,lx_i,ierr)
296: ! Compute initial guess
298: call ApplicationInitialGuess(lx_v(lx_i),ierr)
300: ! Restore vector
302: call VecRestoreArray(X,lx_v,lx_i,ierr)
304: return
305: end
307: ! ---------------------------------------------------------------------
308: !
309: ! ApplicationInitialGuess - Computes initial approximation, called by
310: ! the higher level routine FormInitialGuess().
311: !
312: ! Input Parameter:
313: ! x - local vector data
314: !
315: ! Output Parameters:
316: ! f - local vector data, f(x)
317: ! ierr - error code
318: !
319: ! Notes:
320: ! This routine uses standard Fortran-style computations over a 2-dim array.
321: !
322: subroutine ApplicationInitialGuess(x,ierr)
323: use petscksp
324: implicit none
326: ! Common blocks:
327: PetscReal lambda
328: PetscInt mx,my
329: PetscBool fd_coloring
330: common /params/ lambda,mx,my,fd_coloring
332: ! Input/output variables:
333: PetscScalar x(mx,my)
334: PetscErrorCode ierr
336: ! Local variables:
337: PetscInt i,j
338: PetscReal temp1,temp,hx,hy,one
340: ! Set parameters
342: 0
343: one = 1.0
344: hx = one/(mx-1)
345: hy = one/(my-1)
346: temp1 = lambda/(lambda + one)
348: do 20 j=1,my
349: temp = min(j-1,my-j)*hy
350: do 10 i=1,mx
351: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
352: x(i,j) = 0.0
353: else
354: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
355: endif
356: 10 continue
357: 20 continue
359: return
360: end
362: ! ---------------------------------------------------------------------
363: !
364: ! FormFunction - Evaluates nonlinear function, F(x).
365: !
366: ! Input Parameters:
367: ! snes - the SNES context
368: ! X - input vector
369: ! dummy - optional user-defined context, as set by SNESSetFunction()
370: ! (not used here)
371: !
372: ! Output Parameter:
373: ! F - vector with newly computed function
374: !
375: ! Notes:
376: ! This routine serves as a wrapper for the lower-level routine
377: ! "ApplicationFunction", where the actual computations are
378: ! done using the standard Fortran style of treating the local
379: ! vector data as a multidimensional array over the local mesh.
380: ! This routine merely accesses the local vector data via
381: ! VecGetArray() and VecRestoreArray().
382: !
383: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
384: use petscsnes
385: implicit none
387: ! Input/output variables:
388: SNES snes
389: Vec X,F
390: PetscErrorCode ierr
391: MatFDColoring fdcoloring
393: ! Common blocks:
394: PetscReal lambda
395: PetscInt mx,my
396: PetscBool fd_coloring
397: common /params/ lambda,mx,my,fd_coloring
399: ! Declarations for use with local arrays:
400: PetscScalar lx_v(0:1),lf_v(0:1)
401: PetscOffset lx_i,lf_i
403: PetscInt, pointer :: indices(:)
405: ! Get pointers to vector data.
406: ! - For default PETSc vectors, VecGetArray() returns a pointer to
407: ! the data array. Otherwise, the routine is implementation dependent.
408: ! - You MUST call VecRestoreArray() when you no longer need access to
409: ! the array.
410: ! - Note that the Fortran interface to VecGetArray() differs from the
411: ! C version. See the Fortran chapter of the users manual for details.
413: call VecGetArrayRead(X,lx_v,lx_i,ierr)
414: call VecGetArray(F,lf_v,lf_i,ierr)
416: ! Compute function
418: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)
420: ! Restore vectors
422: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
423: call VecRestoreArray(F,lf_v,lf_i,ierr)
425: call PetscLogFlops(11.0d0*mx*my,ierr)
426: !
427: ! fdcoloring is in the common block and used here ONLY to test the
428: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
429: !
430: if (fd_coloring) then
431: call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
432: print*,'Indices from GetPerturbedColumnsF90'
433: print*,indices
434: call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
435: endif
436: return
437: end
439: ! ---------------------------------------------------------------------
440: !
441: ! ApplicationFunction - Computes nonlinear function, called by
442: ! the higher level routine FormFunction().
443: !
444: ! Input Parameter:
445: ! x - local vector data
446: !
447: ! Output Parameters:
448: ! f - local vector data, f(x)
449: ! ierr - error code
450: !
451: ! Notes:
452: ! This routine uses standard Fortran-style computations over a 2-dim array.
453: !
454: subroutine ApplicationFunction(x,f,ierr)
455: use petscsnes
456: implicit none
458: ! Common blocks:
459: PetscReal lambda
460: PetscInt mx,my
461: PetscBool fd_coloring
462: common /params/ lambda,mx,my,fd_coloring
464: ! Input/output variables:
465: PetscScalar x(mx,my),f(mx,my)
466: PetscErrorCode ierr
468: ! Local variables:
469: PetscScalar two,one,hx,hy
470: PetscScalar hxdhy,hydhx,sc
471: PetscScalar u,uxx,uyy
472: PetscInt i,j
474: 0
475: one = 1.0
476: two = 2.0
477: hx = one/(mx-1)
478: hy = one/(my-1)
479: sc = hx*hy*lambda
480: hxdhy = hx/hy
481: hydhx = hy/hx
483: ! Compute function
485: do 20 j=1,my
486: do 10 i=1,mx
487: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
488: f(i,j) = x(i,j)
489: else
490: u = x(i,j)
491: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
492: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
493: f(i,j) = uxx + uyy - sc*exp(u)
494: endif
495: 10 continue
496: 20 continue
498: return
499: end
501: ! ---------------------------------------------------------------------
502: !
503: ! FormJacobian - Evaluates Jacobian matrix.
504: !
505: ! Input Parameters:
506: ! snes - the SNES context
507: ! x - input vector
508: ! dummy - optional user-defined context, as set by SNESSetJacobian()
509: ! (not used here)
510: !
511: ! Output Parameters:
512: ! jac - Jacobian matrix
513: ! jac_prec - optionally different preconditioning matrix (not used here)
514: ! flag - flag indicating matrix structure
515: !
516: ! Notes:
517: ! This routine serves as a wrapper for the lower-level routine
518: ! "ApplicationJacobian", where the actual computations are
519: ! done using the standard Fortran style of treating the local
520: ! vector data as a multidimensional array over the local mesh.
521: ! This routine merely accesses the local vector data via
522: ! VecGetArray() and VecRestoreArray().
523: !
524: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
525: use petscsnes
526: implicit none
528: ! Input/output variables:
529: SNES snes
530: Vec X
531: Mat jac,jac_prec
532: PetscErrorCode ierr
533: integer dummy
535: ! Common blocks:
536: PetscReal lambda
537: PetscInt mx,my
538: PetscBool fd_coloring
539: common /params/ lambda,mx,my,fd_coloring
541: ! Declarations for use with local array:
542: PetscScalar lx_v(0:1)
543: PetscOffset lx_i
545: ! Get a pointer to vector data
547: call VecGetArrayRead(X,lx_v,lx_i,ierr)
549: ! Compute Jacobian entries
551: call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)
553: ! Restore vector
555: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
557: ! Assemble matrix
559: call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
560: call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
563: return
564: end
566: ! ---------------------------------------------------------------------
567: !
568: ! ApplicationJacobian - Computes Jacobian matrix, called by
569: ! the higher level routine FormJacobian().
570: !
571: ! Input Parameters:
572: ! x - local vector data
573: !
574: ! Output Parameters:
575: ! jac - Jacobian matrix
576: ! jac_prec - optionally different preconditioning matrix (not used here)
577: ! ierr - error code
578: !
579: ! Notes:
580: ! This routine uses standard Fortran-style computations over a 2-dim array.
581: !
582: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
583: use petscsnes
584: implicit none
586: ! Common blocks:
587: PetscReal lambda
588: PetscInt mx,my
589: PetscBool fd_coloring
590: common /params/ lambda,mx,my,fd_coloring
592: ! Input/output variables:
593: PetscScalar x(mx,my)
594: Mat jac,jac_prec
595: PetscErrorCode ierr
597: ! Local variables:
598: PetscInt i,j,row(1),col(5),i1,i5
599: PetscScalar two,one, hx,hy
600: PetscScalar hxdhy,hydhx,sc,v(5)
602: ! Set parameters
604: i1 = 1
605: i5 = 5
606: one = 1.0
607: two = 2.0
608: hx = one/(mx-1)
609: hy = one/(my-1)
610: sc = hx*hy
611: hxdhy = hx/hy
612: hydhx = hy/hx
614: ! Compute entries of the Jacobian matrix
615: ! - Here, we set all entries for a particular row at once.
616: ! - Note that MatSetValues() uses 0-based row and column numbers
617: ! in Fortran as well as in C.
619: do 20 j=1,my
620: row(1) = (j-1)*mx - 1
621: do 10 i=1,mx
622: row(1) = row(1) + 1
623: ! boundary points
624: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
625: call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
626: ! interior grid points
627: else
628: v(1) = -hxdhy
629: v(2) = -hydhx
630: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
631: v(4) = -hydhx
632: v(5) = -hxdhy
633: col(1) = row(1) - mx
634: col(2) = row(1) - 1
635: col(3) = row(1)
636: col(4) = row(1) + 1
637: col(5) = row(1) + mx
638: call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
639: endif
640: 10 continue
641: 20 continue
643: return
644: end
646: !
647: !/*TEST
648: !
649: ! build:
650: ! requires: !single
651: !
652: ! test:
653: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
654: !
655: ! test:
656: ! suffix: 2
657: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
658: !
659: ! test:
660: ! suffix: 3
661: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
662: ! filter: sort -b
663: ! filter_output: sort -b
664: !
665: !TEST*/