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