Actual source code: ex73f90t.F90
petsc-3.11.4 2019-09-28
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*/
26: !
27: ! --------------------------------------------------------------------------
28: !
29: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
30: ! the partial differential equation
31: !
32: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
33: !
34: ! with boundary conditions
35: !
36: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
37: !
38: ! A finite difference approximation with the usual 5-point stencil
39: ! is used to discretize the boundary value problem to obtain a nonlinear
40: ! system of equations.
41: !
42: ! --------------------------------------------------------------------------
43: ! The following define must be used before including any PETSc include files
44: ! into a module or interface. This is because they can't handle declarations
45: ! in them
46: !
47: module petsc_kkt_solver
48: #include <petsc/finclude/petscdm.h>
49: #include <petsc/finclude/petscmat.h>
50: use petscdm
51: use petscmat
52: type petsc_kkt_solver_type
53: DM::da
54: ! temp A block stuff
55: PetscInt mx,my
56: PetscMPIInt rank
57: PetscReal lambda
58: ! Mats
59: Mat::Amat,AmatLin,Bmat,CMat,Dmat
60: IS::isPhi,isLambda
61: end type petsc_kkt_solver_type
63: end module petsc_kkt_solver
65: module petsc_kkt_solver_interfaces
66: use petsc_kkt_solver
68: Interface SNESSetApplicationContext
69: Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
70: #include <petsc/finclude/petscsnes.h>
71: use petscsnes
72: use petsc_kkt_solver
73: SNES:: snesIn
74: type(petsc_kkt_solver_type) ctx
75: PetscErrorCode ierr
76: End Subroutine
77: End Interface SNESSetApplicationContext
79: Interface SNESGetApplicationContext
80: Subroutine SNESGetApplicationContext(snesIn,ctx,ierr)
81: #include <petsc/finclude/petscsnes.h>
82: use petscsnes
83: use petsc_kkt_solver
84: SNES:: snesIn
85: type(petsc_kkt_solver_type), pointer :: ctx
86: PetscErrorCode ierr
87: End Subroutine
88: End Interface SNESGetApplicationContext
89: end module petsc_kkt_solver_interfaces
91: program main
92: #include <petsc/finclude/petscdm.h>
93: #include <petsc/finclude/petscsnes.h>
94: use petscdm
95: use petscdmda
96: use petscsnes
97: use petsc_kkt_solver
98: use petsc_kkt_solver_interfaces
99: implicit none
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Variable declarations
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: !
104: ! Variables:
105: ! mysnes - nonlinear solver
106: ! x, r - solution, residual vectors
107: ! J - Jacobian matrix
108: ! its - iterations for convergence
109: ! Nx, Ny - number of preocessors in x- and y- directions
110: !
111: SNES:: mysnes
112: Vec:: x,r,x2,x1,x1loc,x2loc
113: Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
114: ! Mat:: tmat
115: DM:: daphi,dalam
116: IS:: isglobal(2)
117: PetscErrorCode ierr
118: PetscInt its,N1,N2,i,j,irow,row(1)
119: PetscInt col(1),low,high,lamlow,lamhigh
120: PetscBool flg
121: PetscInt ione,nfour,itwo,nloc,nloclam
122: PetscReal lambda_max,lambda_min
123: type(petsc_kkt_solver_type) solver
124: PetscScalar bval(1),cval(1),one
126: ! Note: Any user-defined Fortran routines (such as FormJacobian)
127: ! MUST be declared as external.
128: external FormInitialGuess,FormJacobian,FormFunction
130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: ! Initialize program
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
134: if (ierr .ne. 0) then
135: print*,'PetscInitialize failed'
136: stop
137: endif
138: call MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr);CHKERRA(ierr)
140: ! Initialize problem parameters
141: lambda_max = 6.81_PETSC_REAL_KIND
142: lambda_min = 0.0
143: solver%lambda = 6.0
144: ione = 1
145: nfour = 4
146: itwo = 2
147: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr)
148: if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then; SETERRA(PETSC_COMM_SELF,1,'Lambda provided with -par is out of range'); 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: call MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
308: call VecDuplicate(x,r,ierr);CHKERRA(ierr)
310: call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr)
312: call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)
314: call SNESSetApplicationContext(mysnes,solver,ierr);CHKERRA(ierr)
316: call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)
318: ! Set function evaluation routine and vector
319: call SNESSetFunction(mysnes,r,FormFunction,solver,ierr);CHKERRA(ierr)
321: call SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr);CHKERRA(ierr)
323: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
324: ! Customize nonlinear solver; set runtime options
325: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
326: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
327: call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr)
329: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
330: ! Evaluate initial guess; then solve nonlinear system.
331: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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(mysnes,x,ierr);CHKERRA(ierr)
338: call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr)
339: call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr)
340: if (solver%rank .eq. 0) then
341: write(6,100) its
342: endif
343: 100 format('Number of SNES iterations = ',i5)
345: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
346: ! Free work space. All PETSc objects should be destroyed when they
347: ! are no longer needed.
348: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: call MatDestroy(KKTmat,ierr);CHKERRA(ierr)
350: call MatDestroy(Amat,ierr);CHKERRA(ierr)
351: call MatDestroy(Dmat,ierr);CHKERRA(ierr)
352: call MatDestroy(Bmat,ierr);CHKERRA(ierr)
353: call MatDestroy(Cmat,ierr);CHKERRA(ierr)
354: call MatDestroy(solver%AmatLin,ierr);CHKERRA(ierr)
355: call ISDestroy(solver%isPhi,ierr);CHKERRA(ierr)
356: call ISDestroy(solver%isLambda,ierr);CHKERRA(ierr)
357: call VecDestroy(x,ierr);CHKERRA(ierr)
358: call VecDestroy(x2,ierr);CHKERRA(ierr)
359: call VecDestroy(x1,ierr);CHKERRA(ierr)
360: call VecDestroy(x1loc,ierr);CHKERRA(ierr)
361: call VecDestroy(x2loc,ierr);CHKERRA(ierr)
362: call VecDestroy(r,ierr);CHKERRA(ierr)
363: call SNESDestroy(mysnes,ierr);CHKERRA(ierr)
364: call DMDestroy(solver%da,ierr);CHKERRA(ierr)
365: call DMDestroy(daphi,ierr);CHKERRA(ierr)
366: call DMDestroy(dalam,ierr);CHKERRA(ierr)
368: call PetscFinalize(ierr)
369: end
371: ! ---------------------------------------------------------------------
372: !
373: ! FormInitialGuess - Forms initial approximation.
374: !
375: ! Input Parameters:
376: ! X - vector
377: !
378: ! Output Parameter:
379: ! X - vector
380: !
381: ! Notes:
382: ! This routine serves as a wrapper for the lower-level routine
383: ! "InitialGuessLocal", where the actual computations are
384: ! done using the standard Fortran style of treating the local
385: ! vector data as a multidimensional array over the local mesh.
386: ! This routine merely handles ghost point scatters and accesses
387: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
388: !
389: subroutine FormInitialGuess(mysnes,Xnest,ierr)
390: #include <petsc/finclude/petscsnes.h>
391: use petscsnes
392: use petsc_kkt_solver
393: use petsc_kkt_solver_interfaces
394: implicit none
395: ! Input/output variables:
396: SNES:: mysnes
397: Vec:: Xnest
398: PetscErrorCode ierr
400: ! Declarations for use with local arrays:
401: type(petsc_kkt_solver_type), pointer:: solver
402: Vec:: Xsub(2)
403: PetscInt:: izero,ione,itwo
405: izero = 0
406: ione = 1
407: itwo = 2
408: 0
409: call SNESGetApplicationContext(mysnes,solver,ierr);CHKERRQ(ierr)
410: call DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
412: call InitialGuessLocal(solver,Xsub(1),ierr);CHKERRQ(ierr)
413: call VecAssemblyBegin(Xsub(1),ierr);CHKERRQ(ierr)
414: call VecAssemblyEnd(Xsub(1),ierr);CHKERRQ(ierr)
416: ! zero out lambda
417: call VecZeroEntries(Xsub(2),ierr);CHKERRQ(ierr)
418: call DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
420: return
421: end subroutine FormInitialGuess
423: ! ---------------------------------------------------------------------
424: !
425: ! InitialGuessLocal - Computes initial approximation, called by
426: ! the higher level routine FormInitialGuess().
427: !
428: ! Input Parameter:
429: ! X1 - local vector data
430: !
431: ! Output Parameters:
432: ! x - local vector data
433: ! ierr - error code
434: !
435: ! Notes:
436: ! This routine uses standard Fortran-style computations over a 2-dim array.
437: !
438: subroutine InitialGuessLocal(solver,X1,ierr)
439: #include <petsc/finclude/petscsys.h>
440: use petscsys
441: use petsc_kkt_solver
442: implicit none
443: ! Input/output variables:
444: type (petsc_kkt_solver_type) solver
445: Vec:: X1
446: PetscErrorCode ierr
448: ! Local variables:
449: PetscInt row,i,j,ione,low,high
450: PetscReal temp1,temp,hx,hy,v
451: PetscReal one
453: ! Set parameters
454: ione = 1
455: 0
456: one = 1.0
457: hx = one/(solver%mx-1)
458: hy = one/(solver%my-1)
459: temp1 = solver%lambda/(solver%lambda + one) + one
461: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
463: do 20 row=low,high-1
464: j = row/solver%mx
465: i = mod(row,solver%mx)
466: temp = min(j,solver%my-j+1)*hy
467: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1 ) then
468: v = 0.0
469: else
470: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
471: endif
472: call VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
473: 20 continue
475: return
476: end subroutine InitialGuessLocal
479: ! ---------------------------------------------------------------------
480: !
481: ! FormJacobian - Evaluates Jacobian matrix.
482: !
483: ! Input Parameters:
484: ! dummy - the SNES context
485: ! x - input vector
486: ! solver - solver data
487: !
488: ! Output Parameters:
489: ! jac - Jacobian matrix
490: ! jac_prec - optionally different preconditioning matrix (not used here)
491: ! flag - flag indicating matrix structure
492: !
493: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
494: #include <petsc/finclude/petscsnes.h>
495: use petscsnes
496: use petsc_kkt_solver
497: implicit none
498: ! Input/output variables:
499: SNES:: dummy
500: Vec:: X
501: Mat:: jac,jac_prec
502: type(petsc_kkt_solver_type) solver
503: PetscErrorCode ierr
505: ! Declarations for use with local arrays:
506: Vec:: Xsub(1)
507: Mat:: Amat
508: PetscInt izero,ione
510: izero = 0
511: ione = 1
513: call DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
515: ! Compute entries for the locally owned part of the Jacobian preconditioner.
516: call MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr);CHKERRQ(ierr)
518: call FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr);CHKERRQ(ierr)
519: call MatDestroy(Amat,ierr);CHKERRQ(ierr) ! discard our reference
520: call DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
522: ! the rest of the matrix is not touched
523: call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
524: call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
525: if (jac .ne. jac_prec) then
526: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
527: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
528: end if
530: ! Tell the matrix we will never add a new nonzero location to the
531: ! matrix. If we do it will generate an error.
532: call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr);CHKERRQ(ierr)
534: return
535: end subroutine FormJacobian
536:
537: ! ---------------------------------------------------------------------
538: !
539: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
540: ! called by the higher level routine FormJacobian().
541: !
542: ! Input Parameters:
543: ! x - local vector data
544: !
545: ! Output Parameters:
546: ! jac - Jacobian preconditioner matrix
547: ! ierr - error code
548: !
549: ! Notes:
550: ! This routine uses standard Fortran-style computations over a 2-dim array.
551: !
552: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
553: #include <petsc/finclude/petscmat.h>
554: use petscmat
555: use petsc_kkt_solver
556: implicit none
557: ! Input/output variables:
558: type (petsc_kkt_solver_type) solver
559: Vec:: X1
560: Mat:: jac
561: logical add_nl_term
562: PetscErrorCode ierr
564: ! Local variables:
565: PetscInt irow,row(1),col(5),i,j
566: PetscInt ione,ifive,low,high,ii
567: PetscScalar two,one,hx,hy,hy2inv
568: PetscScalar hx2inv,sc,v(5)
569: PetscScalar,pointer :: lx_v(:)
571: ! Set parameters
572: ione = 1
573: ifive = 5
574: one = 1.0
575: two = 2.0
576: hx = one/(solver%mx-1)
577: hy = one/(solver%my-1)
578: sc = solver%lambda
579: hx2inv = one/(hx*hx)
580: hy2inv = one/(hy*hy)
582: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
583: call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
585: ii = 0
586: do 20 irow=low,high-1
587: j = irow/solver%mx
588: i = mod(irow,solver%mx)
589: ii = ii + 1 ! one based local index
590: ! boundary points
591: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
592: col(1) = irow
593: row(1) = irow
594: v(1) = one
595: call MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
596: ! interior grid points
597: else
598: v(1) = -hy2inv
599: if(j-1==0) v(1) = 0.0
600: v(2) = -hx2inv
601: if(i-1==0) v(2) = 0.0
602: v(3) = two*(hx2inv + hy2inv)
603: if(add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
604: v(4) = -hx2inv
605: if(i+1==solver%mx-1) v(4) = 0.0
606: v(5) = -hy2inv
607: if(j+1==solver%my-1) v(5) = 0.0
608: col(1) = irow - solver%mx
609: col(2) = irow - 1
610: col(3) = irow
611: col(4) = irow + 1
612: col(5) = irow + solver%mx
613: row(1) = irow
614: call MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
615: endif
616: 20 continue
618: call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
620: return
621: end subroutine FormJacobianLocal
624: ! ---------------------------------------------------------------------
625: !
626: ! FormFunction - Evaluates nonlinear function, F(x).
627: !
628: ! Input Parameters:
629: ! snes - the SNES context
630: ! X - input vector
631: ! dummy - optional user-defined context, as set by SNESSetFunction()
632: ! (not used here)
633: !
634: ! Output Parameter:
635: ! F - function vector
636: !
637: subroutine FormFunction(snesIn,X,F,solver,ierr)
638: #include <petsc/finclude/petscsnes.h>
639: use petscsnes
640: use petsc_kkt_solver
641: implicit none
642: ! Input/output variables:
643: SNES:: snesIn
644: Vec:: X,F
645: PetscErrorCode ierr
646: type (petsc_kkt_solver_type) solver
648: ! Declarations for use with local arrays:
649: Vec:: Xsub(2),Fsub(2)
650: PetscInt izero,ione,itwo
652: ! Scatter ghost points to local vector, using the 2-step process
653: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
654: ! By placing code between these two statements, computations can
655: ! be done while messages are in transition.
656:
657: izero = 0
658: ione = 1
659: itwo = 2
660: call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
661: call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
663: call FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr );CHKERRQ(ierr)
664: call MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
666: ! do rest of operator (linear)
667: call MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr);CHKERRQ(ierr)
668: call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
669: call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr)
671: call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
672: call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
673: return
674: end subroutine formfunction
677: ! ---------------------------------------------------------------------
678: !
679: ! FormFunctionNLTerm - Computes nonlinear function, called by
680: ! the higher level routine FormFunction().
681: !
682: ! Input Parameter:
683: ! x - local vector data
684: !
685: ! Output Parameters:
686: ! f - local vector data, f(x)
687: ! ierr - error code
688: !
689: ! Notes:
690: ! This routine uses standard Fortran-style computations over a 2-dim array.
691: !
692: subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
693: #include <petsc/finclude/petscvec.h>
694: use petscvec
695: use petsc_kkt_solver
696: implicit none
697: ! Input/output variables:
698: type (petsc_kkt_solver_type) solver
699: Vec:: X1,F1
700: PetscErrorCode ierr
701: ! Local variables:
702: PetscScalar one,sc
703: PetscScalar u,v(1)
704: PetscInt i,j,low,high,ii,ione,irow,row(1)
705: PetscScalar,pointer :: lx_v(:)
707: one = 1.0
708: sc = solver%lambda
709: ione = 1
711: call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
712: call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
714: ! Compute function over the locally owned part of the grid
715: ii = 0
716: do 20 irow=low,high-1
717: j = irow/solver%mx
718: i = mod(irow,solver%mx)
719: ii = ii + 1 ! one based local index
720: row(1) = irow
721: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1 ) then
722: v(1) = 0.0
723: else
724: u = lx_v(ii)
725: v(1) = -sc*exp(u)
726: endif
727: call VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
728: 20 continue
730: call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
732: call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr)
733: call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr)
735: 0
736: return
737: end subroutine FormFunctionNLTerm
739: !/*TEST
740: !
741: ! build:
742: ! requires: !single !libpgf90 !complex
743: !
744: ! test:
745: ! nsize: 4
746: ! 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_gamg_est_ksp_type cg
747: !
748: !TEST*/