Actual source code: ex73f90t.F90
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: ! This system (A) is augmented with constraints:
10: !
11: ! A -B * phi = rho
12: ! -C I lam = 0
13: !
14: ! where I is the identity, A is the "normal" Poisson equation, B is the "distributor" of the
15: ! total flux (the first block equation is the flux surface averaging equation). The second
16: ! equation lambda = C * x enforces the surface flux auxiliary equation. B and C have all
17: ! positive entries, areas in C and fraction of area in B.
18: !
20: !
21: ! --------------------------------------------------------------------------
22: !
23: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
24: ! the partial differential equation
25: !
26: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
27: !
28: ! with boundary conditions
29: !
30: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
31: !
32: ! A finite difference approximation with the usual 5-point stencil
33: ! is used to discretize the boundary value problem to obtain a nonlinear
34: ! system of equations.
35: !
36: ! --------------------------------------------------------------------------
37: ! The following define must be used before including any PETSc include files
38: ! into a module or interface. This is because they can't handle declarations
39: ! in them
40: !
41: module petsc_kkt_solver
42: #include <petsc/finclude/petscdm.h>
43: #include <petsc/finclude/petscmat.h>
44: use petscdm
45: use petscmat
46: type petsc_kkt_solver_type
47: DM::da
48: ! temp A block stuff
49: PetscInt mx,my
50: PetscMPIInt rank
51: PetscReal lambda
52: ! Mats
53: Mat::Amat,AmatLin,Bmat,CMat,Dmat
54: IS::isPhi,isLambda
55: end type petsc_kkt_solver_type
57: end module petsc_kkt_solver
59: module petsc_kkt_solver_interfaces
60: use petsc_kkt_solver
62: Interface SNESSetApplicationContext
63: Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
64: #include <petsc/finclude/petscsnes.h>
65: use petscsnes
66: use petsc_kkt_solver
67: SNES:: snesIn
68: type(petsc_kkt_solver_type) ctx
69: PetscErrorCode ierr
70: End Subroutine
71: End Interface SNESSetApplicationContext
73: Interface SNESGetApplicationContext
74: Subroutine SNESGetApplicationContext(snesIn,ctx,ierr)
75: #include <petsc/finclude/petscsnes.h>
76: use petscsnes
77: use petsc_kkt_solver
78: SNES:: snesIn
79: type(petsc_kkt_solver_type), pointer :: ctx
80: PetscErrorCode ierr
81: End Subroutine
82: End Interface SNESGetApplicationContext
83: end module petsc_kkt_solver_interfaces
85: program main
86: #include <petsc/finclude/petscdm.h>
87: #include <petsc/finclude/petscsnes.h>
88: use petscdm
89: use petscdmda
90: use petscsnes
91: use petsc_kkt_solver
92: use petsc_kkt_solver_interfaces
93: implicit none
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Variable declarations
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: !
98: ! Variables:
99: ! mysnes - nonlinear solver
100: ! x, r - solution, residual vectors
101: ! J - Jacobian matrix
102: ! its - iterations for convergence
103: ! Nx, Ny - number of preocessors in x- and y- directions
104: !
105: SNES:: mysnes
106: Vec:: x,r,x2,x1,x1loc,x2loc
107: Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
108: ! Mat:: tmat
109: DM:: daphi,dalam
110: IS:: isglobal(2)
111: PetscErrorCode ierr
112: PetscInt its,N1,N2,i,j,irow,row(1)
113: PetscInt col(1),low,high,lamlow,lamhigh
114: PetscBool flg
115: PetscInt ione,nfour,itwo,nloc,nloclam
116: PetscReal lambda_max,lambda_min
117: type(petsc_kkt_solver_type) solver
118: PetscScalar bval(1),cval(1),one
120: ! Note: Any user-defined Fortran routines (such as FormJacobian)
121: ! MUST be declared as external.
122: external FormInitialGuess,FormJacobian,FormFunction
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Initialize program
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
128: if (ierr .ne. 0) then
129: print*,'PetscInitialize failed'
130: stop
131: endif
132: call MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr);CHKERRA(ierr)
134: ! Initialize problem parameters
135: lambda_max = 6.81_PETSC_REAL_KIND
136: lambda_min = 0.0
137: solver%lambda = 6.0
138: ione = 1
139: nfour = 4
140: itwo = 2
141: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr)
142: if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range'); endif
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: ! Create vector data structures; set function evaluation routine
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: ! just get size
149: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione, &
150: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,daphi,ierr);CHKERRA(ierr)
151: call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)
152: call DMSetUp(daphi,ierr);CHKERRA(ierr)
153: call DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
154: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
155: N1 = solver%my*solver%mx
156: N2 = solver%my
157: flg = .false.
158: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr);CHKERRA(ierr)
159: if (flg) then
160: N2 = 0
161: endif
163: call DMDestroy(daphi,ierr);CHKERRA(ierr)
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: ! Create matrix data structure; set Jacobian evaluation routine
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: call DMShellCreate(PETSC_COMM_WORLD,daphi,ierr);CHKERRA(ierr)
169: call DMSetOptionsPrefix(daphi,'phi_',ierr);CHKERRA(ierr)
170: call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)
172: call VecCreate(PETSC_COMM_WORLD,x1,ierr);CHKERRA(ierr)
173: call VecSetSizes(x1,PETSC_DECIDE,N1,ierr);CHKERRA(ierr)
174: call VecSetFromOptions(x1,ierr);CHKERRA(ierr)
176: call VecGetOwnershipRange(x1,low,high,ierr);CHKERRA(ierr)
177: nloc = high - low
179: call MatCreate(PETSC_COMM_WORLD,Amat,ierr);CHKERRA(ierr)
180: call MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
181: call MatSetUp(Amat,ierr);CHKERRA(ierr)
183: call MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr);CHKERRA(ierr)
184: call MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
185: call MatSetUp(solver%AmatLin,ierr);CHKERRA(ierr)
187: call FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr);CHKERRA(ierr)
188: call MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
189: call MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
191: call DMShellSetGlobalVector(daphi,x1,ierr);CHKERRA(ierr)
192: call DMShellSetMatrix(daphi,Amat,ierr);CHKERRA(ierr)
194: call VecCreate(PETSC_COMM_SELF,x1loc,ierr);CHKERRA(ierr)
195: call VecSetSizes(x1loc,nloc,nloc,ierr);CHKERRA(ierr)
196: call VecSetFromOptions(x1loc,ierr);CHKERRA(ierr)
197: call DMShellSetLocalVector(daphi,x1loc,ierr);CHKERRA(ierr)
199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200: ! Create B, C, & D matrices
201: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: call MatCreate(PETSC_COMM_WORLD,Cmat,ierr);CHKERRA(ierr)
203: call MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr);CHKERRA(ierr)
204: call MatSetUp(Cmat,ierr);CHKERRA(ierr)
205: ! create data for C and B
206: call MatCreate(PETSC_COMM_WORLD,Bmat,ierr);CHKERRA(ierr)
207: call MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr);CHKERRA(ierr)
208: call MatSetUp(Bmat,ierr);CHKERRA(ierr)
209: ! create data for D
210: call MatCreate(PETSC_COMM_WORLD,Dmat,ierr);CHKERRA(ierr)
211: call MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr);CHKERRA(ierr)
212: call MatSetUp(Dmat,ierr);CHKERRA(ierr)
214: call VecCreate(PETSC_COMM_WORLD,x2,ierr);CHKERRA(ierr)
215: call VecSetSizes(x2,PETSC_DECIDE,N2,ierr);CHKERRA(ierr)
216: call VecSetFromOptions(x2,ierr);CHKERRA(ierr)
218: call VecGetOwnershipRange(x2,lamlow,lamhigh,ierr);CHKERRA(ierr)
219: nloclam = lamhigh-lamlow
221: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: ! Set fake B and C
223: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224: one = 1.0
225: if (N2 .gt. 0) then
226: bval(1) = -one/(solver%mx-2)
227: ! cval = -one/(solver%my*solver%mx)
228: cval(1) = -one
229: do 20 irow=low,high-1
230: j = irow/solver%mx ! row in domain
231: i = mod(irow,solver%mx)
232: row(1) = irow
233: col(1) = j
234: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
235: ! no op
236: else
237: call MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr);CHKERRA(ierr)
238: endif
239: row(1) = j
240: call MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
241: 20 continue
242: endif
243: call MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
244: call MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
245: call MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
246: call MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
248: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249: ! Set D (indentity)
250: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: do 30 j=lamlow,lamhigh-1
252: row(1) = j
253: cval(1) = one
254: call MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
255: 30 continue
256: call MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
257: call MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
262: call DMShellCreate(PETSC_COMM_WORLD,dalam,ierr);CHKERRA(ierr)
263: call DMShellSetGlobalVector(dalam,x2,ierr);CHKERRA(ierr)
264: call DMShellSetMatrix(dalam,Dmat,ierr);CHKERRA(ierr)
266: call VecCreate(PETSC_COMM_SELF,x2loc,ierr);CHKERRA(ierr)
267: call VecSetSizes(x2loc,nloclam,nloclam,ierr);CHKERRA(ierr)
268: call VecSetFromOptions(x2loc,ierr);CHKERRA(ierr)
269: call DMShellSetLocalVector(dalam,x2loc,ierr);CHKERRA(ierr)
271: call DMSetOptionsPrefix(dalam,'lambda_',ierr);CHKERRA(ierr)
272: call DMSetFromOptions(dalam,ierr);CHKERRA(ierr)
273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: ! Create field split DA
275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: call DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr);CHKERRA(ierr)
277: call DMCompositeAddDM(solver%da,daphi,ierr);CHKERRA(ierr)
278: call DMCompositeAddDM(solver%da,dalam,ierr);CHKERRA(ierr)
279: call DMSetFromOptions(solver%da,ierr);CHKERRA(ierr)
280: call DMSetUp(solver%da,ierr);CHKERRA(ierr)
281: call DMCompositeGetGlobalISs(solver%da,isglobal,ierr);CHKERRA(ierr)
282: solver%isPhi = isglobal(1)
283: solver%isLambda = isglobal(2)
285: ! cache matrices
286: solver%Amat = Amat
287: solver%Bmat = Bmat
288: solver%Cmat = Cmat
289: solver%Dmat = Dmat
291: matArray(1) = Amat
292: matArray(2) = Bmat
293: matArray(3) = Cmat
294: matArray(4) = Dmat
296: call MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr);CHKERRA(ierr)
297: call MatSetFromOptions(KKTmat,ierr);CHKERRA(ierr)
299: ! Extract global and local vectors from DMDA; then duplicate for remaining
300: ! vectors that are the same types
301: call MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
302: call VecDuplicate(x,r,ierr);CHKERRA(ierr)
304: call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr)
306: call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)
308: call SNESSetApplicationContext(mysnes,solver,ierr);CHKERRA(ierr)
310: call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)
312: ! Set function evaluation routine and vector
313: call SNESSetFunction(mysnes,r,FormFunction,solver,ierr);CHKERRA(ierr)
315: call SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr);CHKERRA(ierr)
317: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318: ! Customize nonlinear solver; set runtime options
319: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
321: call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr)
323: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
324: ! Evaluate initial guess; then solve nonlinear system.
325: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
326: ! Note: The user should initialize the vector, x, with the initial guess
327: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
328: ! to employ an initial guess of zero, the user should explicitly set
329: ! this vector to zero by calling VecSet().
331: call FormInitialGuess(mysnes,x,ierr);CHKERRA(ierr)
332: call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr)
333: call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr)
334: if (solver%rank .eq. 0) then
335: write(6,100) its
336: endif
337: 100 format('Number of SNES iterations = ',i5)
339: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340: ! Free work space. All PETSc objects should be destroyed when they
341: ! are no longer needed.
342: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
343: call MatDestroy(KKTmat,ierr);CHKERRA(ierr)
344: call MatDestroy(Amat,ierr);CHKERRA(ierr)
345: call MatDestroy(Dmat,ierr);CHKERRA(ierr)
346: call MatDestroy(Bmat,ierr);CHKERRA(ierr)
347: call MatDestroy(Cmat,ierr);CHKERRA(ierr)
348: call MatDestroy(solver%AmatLin,ierr);CHKERRA(ierr)
349: call ISDestroy(solver%isPhi,ierr);CHKERRA(ierr)
350: call ISDestroy(solver%isLambda,ierr);CHKERRA(ierr)
351: call VecDestroy(x,ierr);CHKERRA(ierr)
352: call VecDestroy(x2,ierr);CHKERRA(ierr)
353: call VecDestroy(x1,ierr);CHKERRA(ierr)
354: call VecDestroy(x1loc,ierr);CHKERRA(ierr)
355: call VecDestroy(x2loc,ierr);CHKERRA(ierr)
356: call VecDestroy(r,ierr);CHKERRA(ierr)
357: call SNESDestroy(mysnes,ierr);CHKERRA(ierr)
358: call DMDestroy(solver%da,ierr);CHKERRA(ierr)
359: call DMDestroy(daphi,ierr);CHKERRA(ierr)
360: call DMDestroy(dalam,ierr);CHKERRA(ierr)
362: call PetscFinalize(ierr)
363: end
365: ! ---------------------------------------------------------------------
366: !
367: ! FormInitialGuess - Forms initial approximation.
368: !
369: ! Input Parameters:
370: ! X - vector
371: !
372: ! Output Parameter:
373: ! X - vector
374: !
375: ! Notes:
376: ! This routine serves as a wrapper for the lower-level routine
377: ! "InitialGuessLocal", 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 handles ghost point scatters and accesses
381: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
382: !
383: subroutine FormInitialGuess(mysnes,Xnest,ierr)
384: #include <petsc/finclude/petscsnes.h>
385: use petscsnes
386: use petsc_kkt_solver
387: use petsc_kkt_solver_interfaces
388: implicit none
389: ! Input/output variables:
390: SNES:: mysnes
391: Vec:: Xnest
392: PetscErrorCode ierr
394: ! Declarations for use with local arrays:
395: type(petsc_kkt_solver_type), pointer:: solver
396: Vec:: Xsub(2)
397: PetscInt:: izero,ione,itwo
399: izero = 0
400: ione = 1
401: itwo = 2
402: 0
403: call SNESGetApplicationContext(mysnes,solver,ierr);CHKERRQ(ierr)
404: call DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
406: call InitialGuessLocal(solver,Xsub(1),ierr);CHKERRQ(ierr)
407: call VecAssemblyBegin(Xsub(1),ierr);CHKERRQ(ierr)
408: call VecAssemblyEnd(Xsub(1),ierr);CHKERRQ(ierr)
410: ! zero out lambda
411: call VecZeroEntries(Xsub(2),ierr);CHKERRQ(ierr)
412: call DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
414: return
415: end subroutine FormInitialGuess
417: ! ---------------------------------------------------------------------
418: !
419: ! InitialGuessLocal - Computes initial approximation, called by
420: ! the higher level routine FormInitialGuess().
421: !
422: ! Input Parameter:
423: ! X1 - local vector data
424: !
425: ! Output Parameters:
426: ! x - local vector data
427: ! ierr - error code
428: !
429: ! Notes:
430: ! This routine uses standard Fortran-style computations over a 2-dim array.
431: !
432: subroutine InitialGuessLocal(solver,X1,ierr)
433: #include <petsc/finclude/petscsys.h>
434: use petscsys
435: use petsc_kkt_solver
436: implicit none
437: ! Input/output variables:
438: type (petsc_kkt_solver_type) solver
439: Vec:: X1
440: PetscErrorCode ierr
442: ! Local variables:
443: PetscInt row,i,j,ione,low,high
444: PetscReal temp1,temp,hx,hy,v
445: PetscReal one
447: ! Set parameters
448: ione = 1
449: 0
450: one = 1.0
451: hx = one/(solver%mx-1)
452: hy = one/(solver%my-1)
453: temp1 = solver%lambda/(solver%lambda + one) + one
455: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
457: do 20 row=low,high-1
458: j = row/solver%mx
459: i = mod(row,solver%mx)
460: temp = min(j,solver%my-j+1)*hy
461: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
462: v = 0.0
463: else
464: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
465: endif
466: call VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
467: 20 continue
469: return
470: end subroutine InitialGuessLocal
472: ! ---------------------------------------------------------------------
473: !
474: ! FormJacobian - Evaluates Jacobian matrix.
475: !
476: ! Input Parameters:
477: ! dummy - the SNES context
478: ! x - input vector
479: ! solver - solver data
480: !
481: ! Output Parameters:
482: ! jac - Jacobian matrix
483: ! jac_prec - optionally different preconditioning matrix (not used here)
484: ! flag - flag indicating matrix structure
485: !
486: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
487: #include <petsc/finclude/petscsnes.h>
488: use petscsnes
489: use petsc_kkt_solver
490: implicit none
491: ! Input/output variables:
492: SNES:: dummy
493: Vec:: X
494: Mat:: jac,jac_prec
495: type(petsc_kkt_solver_type) solver
496: PetscErrorCode ierr
498: ! Declarations for use with local arrays:
499: Vec:: Xsub(1)
500: Mat:: Amat
501: PetscInt ione
503: ione = 1
505: call DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
507: ! Compute entries for the locally owned part of the Jacobian preconditioner.
508: call MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr);CHKERRQ(ierr)
510: call FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr);CHKERRQ(ierr)
511: call MatDestroy(Amat,ierr);CHKERRQ(ierr) ! discard our reference
512: call DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
514: ! the rest of the matrix is not touched
515: call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
516: call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
517: if (jac .ne. jac_prec) then
518: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
519: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
520: end if
522: ! Tell the matrix we will never add a new nonzero location to the
523: ! matrix. If we do it will generate an error.
524: call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr);CHKERRQ(ierr)
526: return
527: end subroutine FormJacobian
529: ! ---------------------------------------------------------------------
530: !
531: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
532: ! called by the higher level routine FormJacobian().
533: !
534: ! Input Parameters:
535: ! x - local vector data
536: !
537: ! Output Parameters:
538: ! jac - Jacobian preconditioner matrix
539: ! ierr - error code
540: !
541: ! Notes:
542: ! This routine uses standard Fortran-style computations over a 2-dim array.
543: !
544: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
545: #include <petsc/finclude/petscmat.h>
546: use petscmat
547: use petsc_kkt_solver
548: implicit none
549: ! Input/output variables:
550: type (petsc_kkt_solver_type) solver
551: Vec:: X1
552: Mat:: jac
553: logical add_nl_term
554: PetscErrorCode ierr
556: ! Local variables:
557: PetscInt irow,row(1),col(5),i,j
558: PetscInt ione,ifive,low,high,ii
559: PetscScalar two,one,hx,hy,hy2inv
560: PetscScalar hx2inv,sc,v(5)
561: PetscScalar,pointer :: lx_v(:)
563: ! Set parameters
564: ione = 1
565: ifive = 5
566: one = 1.0
567: two = 2.0
568: hx = one/(solver%mx-1)
569: hy = one/(solver%my-1)
570: sc = solver%lambda
571: hx2inv = one/(hx*hx)
572: hy2inv = one/(hy*hy)
574: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
575: call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
577: ii = 0
578: do 20 irow=low,high-1
579: j = irow/solver%mx
580: i = mod(irow,solver%mx)
581: ii = ii + 1 ! one based local index
582: ! boundary points
583: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
584: col(1) = irow
585: row(1) = irow
586: v(1) = one
587: call MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
588: ! interior grid points
589: else
590: v(1) = -hy2inv
591: if (j-1==0) v(1) = 0.0
592: v(2) = -hx2inv
593: if (i-1==0) v(2) = 0.0
594: v(3) = two*(hx2inv + hy2inv)
595: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
596: v(4) = -hx2inv
597: if (i+1==solver%mx-1) v(4) = 0.0
598: v(5) = -hy2inv
599: if (j+1==solver%my-1) v(5) = 0.0
600: col(1) = irow - solver%mx
601: col(2) = irow - 1
602: col(3) = irow
603: col(4) = irow + 1
604: col(5) = irow + solver%mx
605: row(1) = irow
606: call MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
607: endif
608: 20 continue
610: call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
612: return
613: end subroutine FormJacobianLocal
615: ! ---------------------------------------------------------------------
616: !
617: ! FormFunction - Evaluates nonlinear function, F(x).
618: !
619: ! Input Parameters:
620: ! snes - the SNES context
621: ! X - input vector
622: ! dummy - optional user-defined context, as set by SNESSetFunction()
623: ! (not used here)
624: !
625: ! Output Parameter:
626: ! F - function vector
627: !
628: subroutine FormFunction(snesIn,X,F,solver,ierr)
629: #include <petsc/finclude/petscsnes.h>
630: use petscsnes
631: use petsc_kkt_solver
632: implicit none
633: ! Input/output variables:
634: SNES:: snesIn
635: Vec:: X,F
636: PetscErrorCode ierr
637: type (petsc_kkt_solver_type) solver
639: ! Declarations for use with local arrays:
640: Vec:: Xsub(2),Fsub(2)
641: PetscInt itwo
643: ! Scatter ghost points to local vector, using the 2-step process
644: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
645: ! By placing code between these two statements, computations can
646: ! be done while messages are in transition.
648: itwo = 2
649: call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
650: call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
652: call FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr);CHKERRQ(ierr)
653: call MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
655: ! do rest of operator (linear)
656: call MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr);CHKERRQ(ierr)
657: call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
658: call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr)
660: call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
661: call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
662: return
663: end subroutine formfunction
665: ! ---------------------------------------------------------------------
666: !
667: ! FormFunctionNLTerm - Computes nonlinear function, called by
668: ! the higher level routine FormFunction().
669: !
670: ! Input Parameter:
671: ! x - local vector data
672: !
673: ! Output Parameters:
674: ! f - local vector data, f(x)
675: ! ierr - error code
676: !
677: ! Notes:
678: ! This routine uses standard Fortran-style computations over a 2-dim array.
679: !
680: subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
681: #include <petsc/finclude/petscvec.h>
682: use petscvec
683: use petsc_kkt_solver
684: implicit none
685: ! Input/output variables:
686: type (petsc_kkt_solver_type) solver
687: Vec:: X1,F1
688: PetscErrorCode ierr
689: ! Local variables:
690: PetscScalar sc
691: PetscScalar u,v(1)
692: PetscInt i,j,low,high,ii,ione,irow,row(1)
693: PetscScalar,pointer :: lx_v(:)
695: sc = solver%lambda
696: ione = 1
698: call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
699: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
701: ! Compute function over the locally owned part of the grid
702: ii = 0
703: do 20 irow=low,high-1
704: j = irow/solver%mx
705: i = mod(irow,solver%mx)
706: ii = ii + 1 ! one based local index
707: row(1) = irow
708: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
709: v(1) = 0.0
710: else
711: u = lx_v(ii)
712: v(1) = -sc*exp(u)
713: endif
714: call VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
715: 20 continue
717: call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
719: call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr)
720: call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr)
722: 0
723: return
724: end subroutine FormFunctionNLTerm
726: !/*TEST
727: !
728: ! build:
729: ! requires: !single !complex
730: !
731: ! test:
732: ! nsize: 4
733: ! args: -par 5.0 -da_grid_x 10 -da_grid_y 10 -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason -ksp_type fgmres -ksp_norm_type unpreconditioned -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -ksp_monitor_short -fieldsplit_lambda_ksp_type preonly -fieldsplit_lambda_pc_type jacobi -fieldsplit_phi_pc_type gamg -fieldsplit_phi_pc_gamg_esteig_ksp_type cg -fieldsplit_phi_pc_gamg_esteig_ksp_max_it 10 -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0.
734: !
735: !TEST*/