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