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: ! Modified from ex5f90.F by Mike McCourt <mccomic@iit.edu>
10: ! for testing Fortran interface on
11: ! SNESLineSearchSet(), SNESLineSearchSetPreCheck(), SNESLineSearchSetPostCheck()
12: !
13: ! --------------------------------------------------------------------------
14: !
15: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
16: ! the partial differential equation
17: !
18: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
19: !
20: ! with boundary conditions
21: !
22: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
23: !
24: ! A finite difference approximation with the usual 5-point stencil
25: ! is used to discretize the boundary value problem to obtain a nonlinear
26: ! system of equations.
27: !
28: ! The uniprocessor version of this code is snes/examples/tutorials/ex4f.F
29: !
30: ! --------------------------------------------------------------------------
31: ! The following define must be used before including any PETSc include files
32: ! into a module or interface. This is because they can't handle declarations
33: ! in them
34: !
36: module f90module
37: type userctx
38: #include <finclude/petscsysdef.h>
39: #include <finclude/petscvecdef.h>
40: #include <finclude/petscdmdef.h>
41: integer xs,xe,xm,gxs,gxe,gxm
42: integer ys,ye,ym,gys,gye,gym
43: integer mx,my,rank
44: double precision lambda
45: end type userctx
46: contains
47: ! ---------------------------------------------------------------------
48: !
49: ! FormFunction - Evaluates nonlinear function, F(x).
50: !
51: ! Input Parameters:
52: ! snes - the SNES context
53: ! X - input vector
54: ! dummy - optional user-defined context, as set by SNESSetFunction()
55: ! (not used here)
56: !
57: ! Output Parameter:
58: ! F - function vector
59: !
60: ! Notes:
61: ! This routine serves as a wrapper for the lower-level routine
62: ! "FormFunctionLocal", where the actual computations are
63: ! done using the standard Fortran style of treating the local
64: ! vector data as a multidimensional array over the local mesh.
65: ! This routine merely handles ghost point scatters and accesses
66: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
67: !
68: subroutine FormFunction(snes,X,F,user,ierr)
69: implicit none
71: #include <finclude/petscsys.h>
72: #include <finclude/petscvec.h>
73: #include <finclude/petscdmda.h>
74: #include <finclude/petscis.h>
75: #include <finclude/petscmat.h>
76: #include <finclude/petscksp.h>
77: #include <finclude/petscpc.h>
78: #include <finclude/petscsnes.h>
80: #include <finclude/petscvec.h90>
83: ! Input/output variables:
84: SNES snes
85: Vec X,F
86: integer ierr
87: type (userctx) user
88: DM da
90: ! Declarations for use with local arrays:
91: PetscScalar,pointer :: lx_v(:),lf_v(:)
92: Vec localX
94: write(*,*)"Inside FormFunction, user%xm=",user%xm
96: ! Scatter ghost points to local vector, using the 2-step process
97: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
98: ! By placing code between these two statements, computations can
99: ! be done while messages are in transition.
101: call SNESGetDM(snes,da,ierr)
102: call DMGetLocalVector(da,localX,ierr)
103: call DMGlobalToLocalBegin(da,X,INSERT_VALUES, &
104: & localX,ierr)
105: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
107: ! Get a pointer to vector data.
108: ! - For default PETSc vectors, VecGetArray90() returns a pointer to
109: ! the data array. Otherwise, the routine is implementation dependent.
110: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
111: ! the array.
112: ! - Note that the interface to VecGetArrayF90() differs from VecGetArray(),
113: ! and is useable from Fortran-90 Only.
115: call VecGetArrayF90(localX,lx_v,ierr)
116: call VecGetArrayF90(F,lf_v,ierr)
118: ! Compute function over the locally owned part of the grid
120: call FormFunctionLocal(lx_v,lf_v,user,ierr)
122: ! Restore vectors
124: call VecRestoreArrayF90(localX,lx_v,ierr)
125: call VecRestoreArrayF90(F,lf_v,ierr)
127: ! Insert values into global vector
129: call DMRestoreLocalVector(da,localX,ierr)
130: call PetscLogFlops(11.0d0*user%ym*user%xm,ierr)
132: ! call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)
133: ! call VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)
135: return
136: end subroutine formfunction
137: end module f90module
141: program main
142: use f90module
143: implicit none
144: !
145: !
146: #include <finclude/petscsys.h>
147: #include <finclude/petscvec.h>
148: #include <finclude/petscdmda.h>
149: #include <finclude/petscis.h>
150: #include <finclude/petscmat.h>
151: #include <finclude/petscksp.h>
152: #include <finclude/petscpc.h>
153: #include <finclude/petscsnes.h>
154: #include <finclude/petscvec.h90>
155: #include <finclude/petscdmda.h90>
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: ! Variable declarations
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: !
161: ! Variables:
162: ! snes - nonlinear solver
163: ! x, r - solution, residual vectors
164: ! J - Jacobian matrix
165: ! its - iterations for convergence
166: ! Nx, Ny - number of preocessors in x- and y- directions
167: ! matrix_free - flag - 1 indicates matrix-free version
168: !
169: !
170: SNES snes
171: SNESLineSearch linesearch
172: Vec x,r
173: Mat J
174: integer its,matrix_free,flg,ierr
175: double precision lambda_max,lambda_min
176: type (userctx) user
177: PetscBool test_linesearch,test_check
178: DM da
180: ! Note: Any user-defined Fortran routines (such as FormJacobian)
181: ! MUST be declared as external.
183: external FormInitialGuess,FormJacobian,FormLineSearch
184: external FormPostCheck,FormPreCheck
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: ! Initialize program
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
191: call MPI_Comm_rank(PETSC_COMM_WORLD,user%rank,ierr)
193: ! Initialize problem parameters
195: lambda_max = 6.81
196: lambda_min = 0.0
197: user%lambda = 6.0
198: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-par', &
199: & user%lambda,flg,ierr)
200: if (user%lambda .ge. lambda_max .or. user%lambda .le. lambda_min) &
201: & then
202: if (user%rank .eq. 0) write(6,*) 'Lambda is out of range'
203: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
204: endif
207: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: ! Create nonlinear solver context
209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
213: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: ! Create vector data structures; set function evaluation routine
215: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: ! Create distributed array (DMDA) to manage parallel grid and vectors
219: ! This really needs only the star-type stencil, but we use the box
220: ! stencil temporarily.
221: call DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, &
222: & DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX, &
223: & -4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1, &
224: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
225: call DMDAGetInfo(da,PETSC_NULL_INTEGER,user%mx,user%my, &
226: & PETSC_NULL_INTEGER, &
227: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
228: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
229: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
230: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
231: & PETSC_NULL_INTEGER,ierr)
232: 233: !
234: ! Visualize the distribution of the array across the processors
235: !
236: ! call DMView(da,PETSC_VIEWER_DRAW_WORLD,ierr)
238: ! Extract global and local vectors from DMDA; then duplicate for remaining
239: ! vectors that are the same types
240: call SNESSetDM(snes,da,ierr)
241: call DMCreateGlobalVector(da,x,ierr)
242: call VecDuplicate(x,r,ierr)
244: ! Get local grid boundaries (for 2-dimensional DMDA)
246: call DMDAGetCorners(da,user%xs,user%ys,PETSC_NULL_INTEGER, &
247: & user%xm,user%ym,PETSC_NULL_INTEGER,ierr)
248: call DMDAGetGhostCorners(da,user%gxs,user%gys, &
249: & PETSC_NULL_INTEGER,user%gxm,user%gym, &
250: & PETSC_NULL_INTEGER,ierr)
252: ! Here we shift the starting indices up by one so that we can easily
253: ! use the Fortran convention of 1-based indices (rather 0-based indices).
255: user%xs = user%xs+1
256: user%ys = user%ys+1
257: user%gxs = user%gxs+1
258: user%gys = user%gys+1
260: user%ye = user%ys+user%ym-1
261: user%xe = user%xs+user%xm-1
262: user%gye = user%gys+user%gym-1
263: user%gxe = user%gxs+user%gxm-1
265: ! Set function evaluation routine and vector
267: call SNESSetFunction(snes,r,FormFunction,user,ierr)
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: ! Create matrix data structure; set Jacobian evaluation routine
271: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273: ! Set Jacobian matrix data structure and default Jacobian evaluation
274: ! routine. User can override with:
275: ! -snes_fd : default finite differencing approximation of Jacobian
276: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
277: ! (unless user explicitly sets preconditioner)
278: ! -snes_mf_operator : form preconditioning matrix as set by the user,
279: ! but use matrix-free approx for Jacobian-vector
280: ! products within Newton-Krylov method
281: !
282: ! Note: For the parallel case, vectors and matrices MUST be partitioned
283: ! accordingly. When using distributed arrays (DMDAs) to create vectors,
284: ! the DMDAs determine the problem partitioning. We must explicitly
285: ! specify the local matrix dimensions upon its creation for compatibility
286: ! with the vector distribution. Thus, the generic MatCreate() routine
287: ! is NOT sufficient when working with distributed arrays.
288: !
289: ! Note: Here we only approximately preallocate storage space for the
290: ! Jacobian. See the users manual for a discussion of better techniques
291: ! for preallocating matrix memory.
293: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-snes_mf', &
294: & matrix_free,ierr)
295: if (matrix_free .eq. 0) then
296: call DMCreateMatrix(da,MATAIJ,J,ierr)
297: call SNESSetJacobian(snes,J,J,FormJacobian,user,ierr)
298: endif
300: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301: ! Customize nonlinear solver; set runtime options
302: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
306: test_linesearch = PETSC_FALSE307: test_check = PETSC_FALSE308: call PetscOptionsGetBool(PETSC_NULL_CHARACTER,'-test_check',
309: & test_check,PETSC_NULL_INTEGER,ierr)
310: call SNESGetSNESLineSearch(snes, linesearch, ierr)
311: if (test_check.eqv.PETSC_TRUE) then
312: call SNESLineSearchSetPreCheck(linesearch,FormPreCheck,user, &
313: & ierr)
314: call SNESLineSearchSetPostCheck(linesearch,FormPostCheck,user, &
315: & ierr)
316: else
317: call PetscOptionsGetBool(PETSC_NULL_CHARACTER,
318: & '-test_linesearch',test_linesearch,PETSC_NULL_INTEGER,ierr)
319: if (test_linesearch.eqv.PETSC_TRUE) then
320: call SNESLineSearchSetType(linesearch,"shell", ierr)
321: call SNESLineSearchShellSetUserFunc(linesearch, &
322: & FormLineSearch, user, ierr)
323: end if
324: end if
325: call SNESSetFromOptions(snes,ierr)
328: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
329: ! Evaluate initial guess; then solve nonlinear system.
330: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332: ! Note: The user should initialize the vector, x, with the initial guess
333: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
334: ! to employ an initial guess of zero, the user should explicitly set
335: ! this vector to zero by calling VecSet().
337: call FormInitialGuess(snes,x,ierr)
338: write(*,*)"Before SNESSolve"
339: write(*,*)"user%xm=",user%xm
340: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
341: call SNESGetIterationNumber(snes,its,ierr);
342: if (user%rank .eq. 0) then
343: write(6,100) its
344: endif
345: 100 format('Number of SNES iterations = ',i5)
347: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: ! Free work space. All PETSc objects should be destroyed when they
349: ! are no longer needed.
350: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
352: if (matrix_free .eq. 0) call MatDestroy(J,ierr)
353: call VecDestroy(x,ierr)
354: call VecDestroy(r,ierr)
355: call SNESDestroy(snes,ierr)
356: call DMDestroy(da,ierr)
357: call PetscFinalize(ierr)
358: end
360: ! ---------------------------------------------------------------------
361: !
362: ! FormLineSearch - Applies the line search to the step size
363: !
364: subroutine FormLineSearch(linesearch, user, ierr)
366: use f90module
368: #include <finclude/petscsys.h>
369: #include <finclude/petscvec.h>
370: #include <finclude/petscvec.h90>
371: #include <finclude/petscmat.h>
372: #include <finclude/petscmat.h90>
373: #include <finclude/petscksp.h>
374: #include <finclude/petscpc.h>
375: #include <finclude/petscsnes.h>
377: SNES snes
378: type (userctx) user
379: Vec x,f,g,y,w
380: PetscReal fnorm,ynorm,gnorm,xnorm
381: PetscBool flag
382: PetscErrorCode ierr
384: PetscScalar mone
386: write(*,*)"Inside FormLineSearch, user%xm=",user%xm
387: mone = -1.0d0
388: call SNESLineSearchGetSNES(linesearch, snes, ierr)
389: call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
390: call VecNorm(y,NORM_2,ynorm,ierr)
391: call VecAXPY(x,mone,y,ierr)
392: call SNESComputeFunction(snes,x,f,ierr)
393: call VecNorm(f,NORM_2,gnorm,ierr)
394: call VecNorm(x,NORM_2,xnorm,ierr)
395: call SNESLineSearchSetNorms(linesearch,xnorm,gnorm,ynorm,ierr)
396: return
397: end
399: ! ---------------------------------------------------------------------
400: !
401: ! FormInitialGuess - Forms initial approximation.
402: !
403: ! Input Parameters:
404: ! X - vector
405: !
406: ! Output Parameter:
407: ! X - vector
408: !
409: ! Notes:
410: ! This routine serves as a wrapper for the lower-level routine
411: ! "InitialGuessLocal", where the actual computations are
412: ! done using the standard Fortran style of treating the local
413: ! vector data as a multidimensional array over the local mesh.
414: ! This routine merely handles ghost point scatters and accesses
415: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
416: !
417: subroutine FormInitialGuess(snes,X,ierr)
418: use f90module
419: implicit none
421: #include <finclude/petscvec.h90>
422: #include <finclude/petscsys.h>
423: #include <finclude/petscvec.h>
424: #include <finclude/petscdmda.h>
425: #include <finclude/petscis.h>
426: #include <finclude/petscmat.h>
427: #include <finclude/petscksp.h>
428: #include <finclude/petscpc.h>
429: #include <finclude/petscsnes.h>
431: ! Input/output variables:
432: SNES snes
433: Vec X
434: integer ierr
435: DM da
437: ! Declarations for use with local arrays:
438: PetscScalar,pointer :: lx_v(:)
439: Vec localX
440: type (userctx) user
442: 0
444: ! Get a pointer to vector data.
445: ! - For default PETSc vectors, VecGetArray90() returns a pointer to
446: ! the data array. Otherwise, the routine is implementation dependent.
447: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
448: ! the array.
449: ! - Note that the interface to VecGetArrayF90() differs from VecGetArray(),
450: ! and is useable from Fortran-90 Only.
452: call SNESGetApplicationContext(snes,user,ierr)
453: call SNESGetDM(snes,da,ierr)
454: call DMGetLocalVector(da,localX,ierr)
455: call VecGetArrayF90(localX,lx_v,ierr)
457: ! Compute initial guess over the locally owned part of the grid
459: call InitialGuessLocal(user,lx_v,ierr)
461: ! Restore vector
463: call VecRestoreArrayF90(localX,lx_v,ierr)
465: ! Insert values into global vector
467: call DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X,ierr)
468: call DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X,ierr)
469: call DMRestoreLocalVector(da,localX,ierr)
471: return
472: end
474: ! ---------------------------------------------------------------------
475: !
476: ! InitialGuessLocal - Computes initial approximation, called by
477: ! the higher level routine FormInitialGuess().
478: !
479: ! Input Parameter:
480: ! x - local vector data
481: !
482: ! Output Parameters:
483: ! x - local vector data
484: ! ierr - error code
485: !
486: ! Notes:
487: ! This routine uses standard Fortran-style computations over a 2-dim array.
488: !
489: subroutine InitialGuessLocal(user,x,ierr)
490: use f90module
491: implicit none
493: #include <finclude/petscsys.h>
494: #include <finclude/petscvec.h>
495: #include <finclude/petscdmda.h>
496: #include <finclude/petscis.h>
497: #include <finclude/petscmat.h>
498: #include <finclude/petscksp.h>
499: #include <finclude/petscpc.h>
500: #include <finclude/petscsnes.h>
502: ! Input/output variables:
503: type (userctx) user
504: PetscScalar x(user%gxs:user%gxe, &
505: & user%gys:user%gye)
506: integer ierr
508: ! Local variables:
509: integer i,j
510: PetscScalar temp1,temp,hx,hy
511: PetscScalar one
513: ! Set parameters
515: 0
516: one = 1.0
517: hx = one/(dble(user%mx-1))
518: hy = one/(dble(user%my-1))
519: temp1 = user%lambda/(user%lambda + one)
521: do 20 j=user%ys,user%ye
522: temp = dble(min(j-1,user%my-j))*hy
523: do 10 i=user%xs,user%xe
524: if (i .eq. 1 .or. j .eq. 1 &
525: & .or. i .eq. user%mx .or. j .eq. user%my) then
526: x(i,j) = 0.0
527: else
528: x(i,j) = temp1 * &
529: & sqrt(min(dble(min(i-1,user%mx-i)*hx),dble(temp)))
530: endif
531: 10 continue
532: 20 continue
534: return
535: end
537: ! ---------------------------------------------------------------------
538: !
539: ! FormFunctionLocal - Computes nonlinear function, called by
540: ! the higher level routine FormFunction().
541: !
542: ! Input Parameter:
543: ! x - local vector data
544: !
545: ! Output Parameters:
546: ! f - local vector data, f(x)
547: ! ierr - error code
548: !
549: ! Notes:
550: ! This routine uses standard Fortran-style computations over a 2-dim array.
551: !
552: subroutine FormFunctionLocal(x,f,user,ierr)
553: use f90module
555: implicit none
557: ! Input/output variables:
558: type (userctx) user
559: PetscScalar x(user%gxs:user%gxe, &
560: & user%gys:user%gye)
561: PetscScalar f(user%xs:user%xe, &
562: & user%ys:user%ye)
563: integer ierr
565: ! Local variables:
566: PetscScalar two,one,hx,hy,hxdhy,hydhx,sc
567: PetscScalar u,uxx,uyy
568: integer i,j
570: one = 1.0
571: two = 2.0
572: hx = one/dble(user%mx-1)
573: hy = one/dble(user%my-1)
574: sc = hx*hy*user%lambda
575: hxdhy = hx/hy
576: hydhx = hy/hx
578: ! Compute function over the locally owned part of the grid
580: do 20 j=user%ys,user%ye
581: do 10 i=user%xs,user%xe
582: if (i .eq. 1 .or. j .eq. 1 &
583: & .or. i .eq. user%mx .or. j .eq. user%my) then
584: f(i,j) = x(i,j)
585: else
586: u = x(i,j)
587: uxx = hydhx * (two*u &
588: & - x(i-1,j) - x(i+1,j))
589: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
590: f(i,j) = uxx + uyy - sc*exp(u)
591: endif
592: 10 continue
593: 20 continue
595: return
596: end
598: ! ---------------------------------------------------------------------
599: !
600: ! FormJacobian - Evaluates Jacobian matrix.
601: !
602: ! Input Parameters:
603: ! snes - the SNES context
604: ! x - input vector
605: ! dummy - optional user-defined context, as set by SNESSetJacobian()
606: ! (not used here)
607: !
608: ! Output Parameters:
609: ! jac - Jacobian matrix
610: ! jac_prec - optionally different preconditioning matrix (not used here)
611: ! flag - flag indicating matrix structure
612: !
613: ! Notes:
614: ! This routine serves as a wrapper for the lower-level routine
615: ! "FormJacobianLocal", where the actual computations are
616: ! done using the standard Fortran style of treating the local
617: ! vector data as a multidimensional array over the local mesh.
618: ! This routine merely accesses the local vector data via
619: ! VecGetArrayF90() and VecRestoreArrayF90().
620: !
621: ! Notes:
622: ! Due to grid point reordering with DMDAs, we must always work
623: ! with the local grid points, and then transform them to the new
624: ! global numbering with the "ltog" mapping (via DMDAGetGlobalIndicesF90()).
625: ! We cannot work directly with the global numbers for the original
626: ! uniprocessor grid!
627: !
628: ! Two methods are available for imposing this transformation
629: ! when setting matrix entries:
630: ! (A) MatSetValuesLocal(), using the local ordering (including
631: ! ghost points!)
632: ! - Use DMDAGetGlobalIndicesF90() to extract the local-to-global map
633: ! - Associate this map with the matrix by calling
634: ! MatSetLocalToGlobalMapping() once
635: ! - Set matrix entries using the local ordering
636: ! by calling MatSetValuesLocal()
637: ! (B) MatSetValues(), using the global ordering
638: ! - Use DMDAGetGlobalIndicesF90() to extract the local-to-global map
639: ! - Then apply this map explicitly yourself
640: ! - Set matrix entries using the global ordering by calling
641: ! MatSetValues()
642: ! Option (A) seems cleaner/easier in many cases, and is the procedure
643: ! used in this example.
644: !
645: subroutine FormJacobian(snes,X,jac,jac_prec,flag,user,ierr)
646: use f90module
647: implicit none
649: #include <finclude/petscsys.h>
650: #include <finclude/petscvec.h>
651: #include <finclude/petscdmda.h>
652: #include <finclude/petscis.h>
653: #include <finclude/petscmat.h>
654: #include <finclude/petscksp.h>
655: #include <finclude/petscpc.h>
656: #include <finclude/petscsnes.h>
658: #include <finclude/petscvec.h90>
660: ! Input/output variables:
661: SNES snes
662: Vec X
663: Mat jac,jac_prec
664: MatStructure flag
665: type(userctx) user
666: integer ierr
667: DM da
669: ! Declarations for use with local arrays:
670: PetscScalar,pointer :: lx_v(:)
671: Vec localX
673: ! Scatter ghost points to local vector, using the 2-step process
674: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd()
675: ! Computations can be done while messages are in transition,
676: ! by placing code between these two statements.
678: call SNESGetDM(snes,da,ierr)
679: call DMGetLocalVector(da,localX,ierr)
680: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX, &
681: & ierr)
682: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
684: ! Get a pointer to vector data
686: call VecGetArrayF90(localX,lx_v,ierr)
688: ! Compute entries for the locally owned part of the Jacobian.
690: call FormJacobianLocal(lx_v,jac,jac_prec,user,ierr)
692: ! Assemble matrix, using the 2-step process:
693: ! MatAssemblyBegin(), MatAssemblyEnd()
694: ! Computations can be done while messages are in transition,
695: ! by placing code between these two statements.
697: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
698: call VecRestoreArrayF90(localX,lx_v,ierr)
699: call DMRestoreLocalVector(da,localX,ierr)
700: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
702: ! Set flag to indicate that the Jacobian matrix retains an identical
703: ! nonzero structure throughout all nonlinear iterations (although the
704: ! values of the entries change). Thus, we can save some work in setting
705: ! up the preconditioner (e.g., no need to redo symbolic factorization for
706: ! ILU/ICC preconditioners).
707: ! - If the nonzero structure of the matrix is different during
708: ! successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
709: ! must be used instead. If you are unsure whether the matrix
710: ! structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
711: ! - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
712: ! believes your assertion and does not check the structure
713: ! of the matrix. If you erroneously claim that the structure
714: ! is the same when it actually is not, the new preconditioner
715: ! will not function correctly. Thus, use this optimization
716: ! feature with caution!
718: flag = SAME_NONZERO_PATTERN
720: ! Tell the matrix we will never add a new nonzero location to the
721: ! matrix. If we do it will generate an error.
723: ! call MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,ierr)
725: return
726: end
728: ! ---------------------------------------------------------------------
729: !
730: ! FormJacobianLocal - Computes Jacobian matrix, called by
731: ! the higher level routine FormJacobian().
732: !
733: ! Input Parameters:
734: ! x - local vector data
735: !
736: ! Output Parameters:
737: ! jac - Jacobian matrix
738: ! jac_prec - optionally different preconditioning matrix (not used here)
739: ! ierr - error code
740: !
741: ! Notes:
742: ! This routine uses standard Fortran-style computations over a 2-dim array.
743: !
744: ! Notes:
745: ! Due to grid point reordering with DMDAs, we must always work
746: ! with the local grid points, and then transform them to the new
747: ! global numbering with the "ltog" mapping (via DMDAGetGlobalIndicesF90()).
748: ! We cannot work directly with the global numbers for the original
749: ! uniprocessor grid!
750: !
751: ! Two methods are available for imposing this transformation
752: ! when setting matrix entries:
753: ! (A) MatSetValuesLocal(), using the local ordering (including
754: ! ghost points!)
755: ! - Use DMDAGetGlobalIndicesF90() to extract the local-to-global map
756: ! - Associate this map with the matrix by calling
757: ! MatSetLocalToGlobalMapping() once
758: ! - Set matrix entries using the local ordering
759: ! by calling MatSetValuesLocal()
760: ! (B) MatSetValues(), using the global ordering
761: ! - Use DMDAGetGlobalIndicesF90() to extract the local-to-global map
762: ! - Then apply this map explicitly yourself
763: ! - Set matrix entries using the global ordering by calling
764: ! MatSetValues()
765: ! Option (A) seems cleaner/easier in many cases, and is the procedure
766: ! used in this example.
767: !
768: subroutine FormJacobianLocal(x,jac,jac_prec,user,ierr)
769: use f90module
770: implicit none
772: #include <finclude/petscsys.h>
773: #include <finclude/petscvec.h>
774: #include <finclude/petscdmda.h>
775: #include <finclude/petscis.h>
776: #include <finclude/petscmat.h>
777: #include <finclude/petscksp.h>
778: #include <finclude/petscpc.h>
779: #include <finclude/petscsnes.h>
781: ! Input/output variables:
782: type (userctx) user
783: PetscScalar x(user%gxs:user%gxe, &
784: & user%gys:user%gye)
785: Mat jac,jac_prec
786: integer ierr
788: ! Local variables:
789: integer row,col(5),i,j
790: PetscScalar two,one,hx,hy,hxdhy
791: PetscScalar hydhx,sc,v(5)
793: ! Set parameters
795: one = 1.0
796: two = 2.0
797: hx = one/dble(user%mx-1)
798: hy = one/dble(user%my-1)
799: sc = hx*hy
800: hxdhy = hx/hy
801: hydhx = hy/hx
804: ! Compute entries for the locally owned part of the Jacobian.
805: ! - Currently, all PETSc parallel matrix formats are partitioned by
806: ! contiguous chunks of rows across the processors.
807: ! - Each processor needs to insert only elements that it owns
808: ! locally (but any non-local elements will be sent to the
809: ! appropriate processor during matrix assembly).
810: ! - Here, we set all entries for a particular row at once.
811: ! - We can set matrix entries either using either
812: ! MatSetValuesLocal() or MatSetValues(), as discussed above.
813: ! - Note that MatSetValues() uses 0-based row and column numbers
814: ! in Fortran as well as in C.
816: do 20 j=user%ys,user%ye
817: row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1
818: do 10 i=user%xs,user%xe
819: row = row + 1
820: ! boundary points
821: if (i .eq. 1 .or. j .eq. 1 &
822: & .or. i .eq. user%mx .or. j .eq. user%my) then
823: col(1) = row
824: v(1) = one
825: call MatSetValuesLocal(jac,1,row,1,col,v, &
826: & INSERT_VALUES,ierr)
827: ! interior grid points
828: else
829: v(1) = -hxdhy
830: v(2) = -hydhx
831: v(3) = two*(hydhx + hxdhy) &
832: & - sc*user%lambda*exp(x(i,j))
833: v(4) = -hydhx
834: v(5) = -hxdhy
835: col(1) = row - user%gxm
836: col(2) = row - 1
837: col(3) = row
838: col(4) = row + 1
839: col(5) = row + user%gxm
840: call MatSetValuesLocal(jac,1,row,5,col,v, &
841: & INSERT_VALUES,ierr)
842: endif
843: 10 continue
844: 20 continue
846: return
847: end
849: subroutine FormPreCheck(snes,X,Y,changed_Y,user,ierr)
850: use f90module
852: SNES snes
853: Vec X,Y
854: type (userctx) user
855: PetscBool changed_Y
856: PetscErrorCode ierr
858: write(*,*)"Inside formPreCheck, user%xm=",user%xm
859: end subroutine formPreCheck
861: subroutine FormPostCheck(snes,X,Y,W,changed_Y,changed_W,user,ierr)
862: use f90module
864: SNES snes
865: Vec X,Y,W
866: type (userctx) user
867: PetscBool changed_Y,changed_W
868: PetscErrorCode ierr
870: write(*,*)"Inside formPostCheck, user%xm=",user%xm
871: end subroutine formPostCheck