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 ex73f90tmodule
42: #include <petsc/finclude/petscdm.h>
43: #include <petsc/finclude/petscmat.h>
44: use petscdm
45: use petscmat
46: type ex73f90tmodule_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 ex73f90tmodule_type
57: end module ex73f90tmodule
59: module ex73f90tmodule_interfaces
60: use ex73f90tmodule
62: Interface SNESSetApplicationContext
63: Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
64: #include <petsc/finclude/petscsnes.h>
65: use petscsnes
66: use ex73f90tmodule
67: SNES:: snesIn
68: type(ex73f90tmodule_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 ex73f90tmodule
78: SNES:: snesIn
79: type(ex73f90tmodule_type), pointer :: ctx
80: PetscErrorCode ierr
81: End Subroutine
82: End Interface SNESGetApplicationContext
83: end module ex73f90tmodule_interfaces
85: subroutine MyObjective(snes, x, result, ctx, ierr )
86: #include <petsc/finclude/petsc.h>
87: use petsc
88: implicit none
89: PetscInt ctx
90: Vec x, f
91: SNES snes
92: PetscErrorCode ierr
93: PetscScalar result
94: PetscReal fnorm
96: PetscCall(VecDuplicate(x,f,ierr))
97: PetscCall(SNESComputeFunction(snes,x,f,ierr))
98: PetscCall(VecNorm(f,NORM_2,fnorm,ierr))
99: result = .5*fnorm*fnorm
100: PetscCall(VecDestroy(f,ierr))
101: end subroutine MyObjective
103: program main
104: #include <petsc/finclude/petscdm.h>
105: #include <petsc/finclude/petscsnes.h>
106: use petscdm
107: use petscdmda
108: use petscsnes
109: use ex73f90tmodule
110: use ex73f90tmodule_interfaces
111: implicit none
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: ! Variable declarations
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: !
116: ! Variables:
117: ! mysnes - nonlinear solver
118: ! x, r - solution, residual vectors
119: ! J - Jacobian matrix
120: ! its - iterations for convergence
121: ! Nx, Ny - number of preocessors in x- and y- directions
122: !
123: SNES:: mysnes
124: Vec:: x,r,x2,x1,x1loc,x2loc
125: Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
126: ! Mat:: tmat
127: DM:: daphi,dalam
128: IS:: isglobal(2)
129: PetscErrorCode ierr
130: PetscInt its,N1,N2,i,j,irow,row(1)
131: PetscInt col(1),low,high,lamlow,lamhigh
132: PetscBool flg
133: PetscInt ione,nfour,itwo,nloc,nloclam
134: PetscReal lambda_max,lambda_min
135: type(ex73f90tmodule_type) solver
136: PetscScalar bval(1),cval(1),one
137: PetscBool useobjective
139: ! Note: Any user-defined Fortran routines (such as FormJacobian)
140: ! MUST be declared as external.
141: external FormInitialGuess,FormJacobian,FormFunction,MyObjective
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Initialize program
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: PetscCallA(PetscInitialize(ierr))
147: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr))
149: ! Initialize problem parameters
150: lambda_max = 6.81_PETSC_REAL_KIND
151: lambda_min = 0.0
152: solver%lambda = 6.0
153: ione = 1
154: nfour = 4
155: itwo = 2
156: useobjective = PETSC_FALSE
157: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr))
158: PetscCheckA(solver%lambda .le. lambda_max .and. solver%lambda .ge. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range')
159: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-objective', useobjective,flg,ierr))
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! Create vector data structures; set function evaluation routine
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: ! just get size
166: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,daphi,ierr))
167: PetscCallA(DMSetFromOptions(daphi,ierr))
168: PetscCallA(DMSetUp(daphi,ierr))
169: PetscCallA(DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,ierr))
170: N1 = solver%my*solver%mx
171: N2 = solver%my
172: flg = .false.
173: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr))
174: if (flg) then
175: N2 = 0
176: endif
178: PetscCallA(DMDestroy(daphi,ierr))
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: ! Create matrix data structure; set Jacobian evaluation routine
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,daphi,ierr))
184: PetscCallA(DMSetOptionsPrefix(daphi,'phi_',ierr))
185: PetscCallA(DMSetFromOptions(daphi,ierr))
187: PetscCallA(VecCreate(PETSC_COMM_WORLD,x1,ierr))
188: PetscCallA(VecSetSizes(x1,PETSC_DECIDE,N1,ierr))
189: PetscCallA(VecSetFromOptions(x1,ierr))
191: PetscCallA(VecGetOwnershipRange(x1,low,high,ierr))
192: nloc = high - low
194: PetscCallA(MatCreate(PETSC_COMM_WORLD,Amat,ierr))
195: PetscCallA(MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
196: PetscCallA(MatSetUp(Amat,ierr))
198: PetscCallA(MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr))
199: PetscCallA(MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr))
200: PetscCallA(MatSetUp(solver%AmatLin,ierr))
202: PetscCallA(FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr))
203: PetscCallA(MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
204: PetscCallA(MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr))
206: PetscCallA(DMShellSetGlobalVector(daphi,x1,ierr))
207: PetscCallA(DMShellSetMatrix(daphi,Amat,ierr))
209: PetscCallA(VecCreate(PETSC_COMM_SELF,x1loc,ierr))
210: PetscCallA(VecSetSizes(x1loc,nloc,nloc,ierr))
211: PetscCallA(VecSetFromOptions(x1loc,ierr))
212: PetscCallA(DMShellSetLocalVector(daphi,x1loc,ierr))
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: ! Create B, C, & D matrices
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: PetscCallA(MatCreate(PETSC_COMM_WORLD,Cmat,ierr))
218: PetscCallA(MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr))
219: PetscCallA(MatSetUp(Cmat,ierr))
220: ! create data for C and B
221: PetscCallA(MatCreate(PETSC_COMM_WORLD,Bmat,ierr))
222: PetscCallA(MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr))
223: PetscCallA(MatSetUp(Bmat,ierr))
224: ! create data for D
225: PetscCallA(MatCreate(PETSC_COMM_WORLD,Dmat,ierr))
226: PetscCallA(MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr))
227: PetscCallA(MatSetUp(Dmat,ierr))
229: PetscCallA(VecCreate(PETSC_COMM_WORLD,x2,ierr))
230: PetscCallA(VecSetSizes(x2,PETSC_DECIDE,N2,ierr))
231: PetscCallA(VecSetFromOptions(x2,ierr))
233: PetscCallA(VecGetOwnershipRange(x2,lamlow,lamhigh,ierr))
234: nloclam = lamhigh-lamlow
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: ! Set fake B and C
238: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: one = 1.0
240: if (N2 .gt. 0) then
241: bval(1) = -one/(solver%mx-2)
242: ! cval = -one/(solver%my*solver%mx)
243: cval(1) = -one
244: do 20 irow=low,high-1
245: j = irow/solver%mx ! row in domain
246: i = mod(irow,solver%mx)
247: row(1) = irow
248: col(1) = j
249: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
250: ! no op
251: else
252: PetscCallA(MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr))
253: endif
254: row(1) = j
255: PetscCallA(MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
256: 20 continue
257: endif
258: PetscCallA(MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr))
259: PetscCallA(MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr))
260: PetscCallA(MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr))
261: PetscCallA(MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr))
263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: ! Set D (identity)
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266: do 30 j=lamlow,lamhigh-1
267: row(1) = j
268: cval(1) = one
269: PetscCallA(MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr))
270: 30 continue
271: PetscCallA(MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr))
272: PetscCallA(MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr))
274: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: ! DM for lambda (dalam) : temp driver for A block, setup A block solver data
276: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: PetscCallA(DMShellCreate(PETSC_COMM_WORLD,dalam,ierr))
278: PetscCallA(DMShellSetGlobalVector(dalam,x2,ierr))
279: PetscCallA(DMShellSetMatrix(dalam,Dmat,ierr))
281: PetscCallA(VecCreate(PETSC_COMM_SELF,x2loc,ierr))
282: PetscCallA(VecSetSizes(x2loc,nloclam,nloclam,ierr))
283: PetscCallA(VecSetFromOptions(x2loc,ierr))
284: PetscCallA(DMShellSetLocalVector(dalam,x2loc,ierr))
286: PetscCallA(DMSetOptionsPrefix(dalam,'lambda_',ierr))
287: PetscCallA(DMSetFromOptions(dalam,ierr))
288: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: ! Create field split DA
290: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291: PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr))
292: PetscCallA(DMCompositeAddDM(solver%da,daphi,ierr))
293: PetscCallA(DMCompositeAddDM(solver%da,dalam,ierr))
294: PetscCallA(DMSetFromOptions(solver%da,ierr))
295: PetscCallA(DMSetUp(solver%da,ierr))
296: PetscCallA(DMCompositeGetGlobalISs(solver%da,isglobal,ierr))
297: solver%isPhi = isglobal(1)
298: solver%isLambda = isglobal(2)
300: ! cache matrices
301: solver%Amat = Amat
302: solver%Bmat = Bmat
303: solver%Cmat = Cmat
304: solver%Dmat = Dmat
306: matArray(1) = Amat
307: matArray(2) = Bmat
308: matArray(3) = Cmat
309: matArray(4) = Dmat
311: PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr))
312: PetscCallA(MatSetFromOptions(KKTmat,ierr))
314: ! Extract global and local vectors from DMDA; then duplicate for remaining
315: ! vectors that are the same types
316: PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr))
317: PetscCallA(VecDuplicate(x,r,ierr))
319: PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr))
321: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
323: PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr))
325: PetscCallA(SNESSetDM(mysnes,solver%da,ierr))
327: ! Set function evaluation routine and vector
328: PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver,ierr))
329: if (useobjective .eqv. PETSC_TRUE) then
330: PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
331: endif
332: PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr))
334: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335: ! Customize nonlinear solver; set runtime options
336: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
338: PetscCallA(SNESSetFromOptions(mysnes,ierr))
340: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
341: ! Evaluate initial guess; then solve nonlinear system.
342: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
343: ! Note: The user should initialize the vector, x, with the initial guess
344: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
345: ! to employ an initial guess of zero, the user should explicitly set
346: ! this vector to zero by calling VecSet().
348: PetscCallA(FormInitialGuess(mysnes,x,ierr))
349: PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr))
350: PetscCallA(SNESGetIterationNumber(mysnes,its,ierr))
351: if (solver%rank .eq. 0) then
352: write(6,100) its
353: endif
354: 100 format('Number of SNES iterations = ',i5)
356: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: ! Free work space. All PETSc objects should be destroyed when they
358: ! are no longer needed.
359: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
360: PetscCallA(MatDestroy(KKTmat,ierr))
361: PetscCallA(MatDestroy(Amat,ierr))
362: PetscCallA(MatDestroy(Dmat,ierr))
363: PetscCallA(MatDestroy(Bmat,ierr))
364: PetscCallA(MatDestroy(Cmat,ierr))
365: PetscCallA(MatDestroy(solver%AmatLin,ierr))
366: PetscCallA(ISDestroy(solver%isPhi,ierr))
367: PetscCallA(ISDestroy(solver%isLambda,ierr))
368: PetscCallA(VecDestroy(x,ierr))
369: PetscCallA(VecDestroy(x2,ierr))
370: PetscCallA(VecDestroy(x1,ierr))
371: PetscCallA(VecDestroy(x1loc,ierr))
372: PetscCallA(VecDestroy(x2loc,ierr))
373: PetscCallA(VecDestroy(r,ierr))
374: PetscCallA(SNESDestroy(mysnes,ierr))
375: PetscCallA(DMDestroy(solver%da,ierr))
376: PetscCallA(DMDestroy(daphi,ierr))
377: PetscCallA(DMDestroy(dalam,ierr))
379: PetscCallA(PetscFinalize(ierr))
380: end
382: ! ---------------------------------------------------------------------
383: !
384: ! FormInitialGuess - Forms initial approximation.
385: !
386: ! Input Parameters:
387: ! X - vector
388: !
389: ! Output Parameter:
390: ! X - vector
391: !
392: ! Notes:
393: ! This routine serves as a wrapper for the lower-level routine
394: ! "InitialGuessLocal", where the actual computations are
395: ! done using the standard Fortran style of treating the local
396: ! vector data as a multidimensional array over the local mesh.
397: ! This routine merely handles ghost point scatters and accesses
398: ! the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
399: !
400: subroutine FormInitialGuess(mysnes,Xnest,ierr)
401: #include <petsc/finclude/petscsnes.h>
402: use petscsnes
403: use ex73f90tmodule
404: use ex73f90tmodule_interfaces
405: implicit none
406: ! Input/output variables:
407: SNES:: mysnes
408: Vec:: Xnest
409: PetscErrorCode ierr
411: ! Declarations for use with local arrays:
412: type(ex73f90tmodule_type), pointer:: solver
413: Vec:: Xsub(2)
414: PetscInt:: izero,ione,itwo
416: izero = 0
417: ione = 1
418: itwo = 2
419: ierr = 0
420: PetscCall(SNESGetApplicationContext(mysnes,solver,ierr))
421: PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
423: PetscCall(InitialGuessLocal(solver,Xsub(1),ierr))
424: PetscCall(VecAssemblyBegin(Xsub(1),ierr))
425: PetscCall(VecAssemblyEnd(Xsub(1),ierr))
427: ! zero out lambda
428: PetscCall(VecZeroEntries(Xsub(2),ierr))
429: PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
431: end subroutine FormInitialGuess
433: ! ---------------------------------------------------------------------
434: !
435: ! InitialGuessLocal - Computes initial approximation, called by
436: ! the higher level routine FormInitialGuess().
437: !
438: ! Input Parameter:
439: ! X1 - local vector data
440: !
441: ! Output Parameters:
442: ! x - local vector data
443: ! ierr - error code
444: !
445: ! Notes:
446: ! This routine uses standard Fortran-style computations over a 2-dim array.
447: !
448: subroutine InitialGuessLocal(solver,X1,ierr)
449: #include <petsc/finclude/petscsys.h>
450: use petscsys
451: use ex73f90tmodule
452: implicit none
453: ! Input/output variables:
454: type (ex73f90tmodule_type) solver
455: Vec:: X1
456: PetscErrorCode ierr
458: ! Local variables:
459: PetscInt row,i,j,ione,low,high
460: PetscReal temp1,temp,hx,hy,v
461: PetscReal one
463: ! Set parameters
464: ione = 1
465: ierr = 0
466: one = 1.0
467: hx = one/(solver%mx-1)
468: hy = one/(solver%my-1)
469: temp1 = solver%lambda/(solver%lambda + one) + one
471: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
473: do 20 row=low,high-1
474: j = row/solver%mx
475: i = mod(row,solver%mx)
476: temp = min(j,solver%my-j+1)*hy
477: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
478: v = 0.0
479: else
480: v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
481: endif
482: PetscCall(VecSetValues(X1,ione,[row],[v],INSERT_VALUES,ierr))
483: 20 continue
485: end subroutine InitialGuessLocal
487: ! ---------------------------------------------------------------------
488: !
489: ! FormJacobian - Evaluates Jacobian matrix.
490: !
491: ! Input Parameters:
492: ! dummy - the SNES context
493: ! x - input vector
494: ! solver - solver data
495: !
496: ! Output Parameters:
497: ! jac - Jacobian matrix
498: ! jac_prec - optionally different preconditioning matrix (not used here)
499: ! flag - flag indicating matrix structure
500: !
501: subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
502: #include <petsc/finclude/petscsnes.h>
503: use petscsnes
504: use ex73f90tmodule
505: implicit none
506: ! Input/output variables:
507: SNES:: dummy
508: Vec:: X
509: Mat:: jac,jac_prec
510: type(ex73f90tmodule_type) solver
511: PetscErrorCode ierr
513: ! Declarations for use with local arrays:
514: Vec:: Xsub(1)
515: Mat:: Amat
516: PetscInt ione
518: ione = 1
520: PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
522: ! Compute entries for the locally owned part of the Jacobian preconditioner.
523: PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr))
525: PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr))
526: PetscCall(MatDestroy(Amat,ierr)) ! discard our reference
527: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
529: ! the rest of the matrix is not touched
530: PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
531: PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
532: if (jac .ne. jac_prec) then
533: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
534: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
535: end if
537: ! Tell the matrix we will never add a new nonzero location to the
538: ! matrix. If we do it will generate an error.
539: PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
541: end subroutine FormJacobian
543: ! ---------------------------------------------------------------------
544: !
545: ! FormJacobianLocal - Computes Jacobian preconditioner matrix,
546: ! called by the higher level routine FormJacobian().
547: !
548: ! Input Parameters:
549: ! x - local vector data
550: !
551: ! Output Parameters:
552: ! jac - Jacobian preconditioner matrix
553: ! ierr - error code
554: !
555: ! Notes:
556: ! This routine uses standard Fortran-style computations over a 2-dim array.
557: !
558: subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
559: #include <petsc/finclude/petscmat.h>
560: use petscmat
561: use ex73f90tmodule
562: implicit none
563: ! Input/output variables:
564: type (ex73f90tmodule_type) solver
565: Vec:: X1
566: Mat:: jac
567: logical add_nl_term
568: PetscErrorCode ierr
570: ! Local variables:
571: PetscInt irow,row(1),col(5),i,j
572: PetscInt ione,ifive,low,high,ii
573: PetscScalar two,one,hx,hy,hy2inv
574: PetscScalar hx2inv,sc,v(5)
575: PetscScalar,pointer :: lx_v(:)
577: ! Set parameters
578: ione = 1
579: ifive = 5
580: one = 1.0
581: two = 2.0
582: hx = one/(solver%mx-1)
583: hy = one/(solver%my-1)
584: sc = solver%lambda
585: hx2inv = one/(hx*hx)
586: hy2inv = one/(hy*hy)
588: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
589: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
591: ii = 0
592: do 20 irow=low,high-1
593: j = irow/solver%mx
594: i = mod(irow,solver%mx)
595: ii = ii + 1 ! one based local index
596: ! boundary points
597: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
598: col(1) = irow
599: row(1) = irow
600: v(1) = one
601: PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr))
602: ! interior grid points
603: else
604: v(1) = -hy2inv
605: if (j-1==0) v(1) = 0.0
606: v(2) = -hx2inv
607: if (i-1==0) v(2) = 0.0
608: v(3) = two*(hx2inv + hy2inv)
609: if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
610: v(4) = -hx2inv
611: if (i+1==solver%mx-1) v(4) = 0.0
612: v(5) = -hy2inv
613: if (j+1==solver%my-1) v(5) = 0.0
614: col(1) = irow - solver%mx
615: col(2) = irow - 1
616: col(3) = irow
617: col(4) = irow + 1
618: col(5) = irow + solver%mx
619: row(1) = irow
620: PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr))
621: endif
622: 20 continue
624: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
626: end subroutine FormJacobianLocal
628: ! ---------------------------------------------------------------------
629: !
630: ! FormFunction - Evaluates nonlinear function, F(x).
631: !
632: ! Input Parameters:
633: ! snes - the SNES context
634: ! X - input vector
635: ! dummy - optional user-defined context, as set by SNESSetFunction()
636: ! (not used here)
637: !
638: ! Output Parameter:
639: ! F - function vector
640: !
641: subroutine FormFunction(snesIn,X,F,solver,ierr)
642: #include <petsc/finclude/petscsnes.h>
643: use petscsnes
644: use ex73f90tmodule
645: implicit none
646: ! Input/output variables:
647: SNES:: snesIn
648: Vec:: X,F
649: PetscErrorCode ierr
650: type (ex73f90tmodule_type) solver
652: ! Declarations for use with local arrays:
653: Vec:: Xsub(2),Fsub(2)
654: PetscInt itwo
656: ! Scatter ghost points to local vector, using the 2-step process
657: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
658: ! By placing code between these two statements, computations can
659: ! be done while messages are in transition.
661: itwo = 2
662: PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
663: PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER_ARRAY,Fsub,ierr))
665: PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr))
666: PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
668: ! do rest of operator (linear)
669: PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr))
670: PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
671: PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
673: PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr))
674: PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER_ARRAY,Fsub,ierr))
675: 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 ex73f90tmodule
696: implicit none
697: ! Input/output variables:
698: type (ex73f90tmodule_type) solver
699: Vec:: X1,F1
700: PetscErrorCode ierr
701: ! Local variables:
702: PetscScalar sc
703: PetscScalar u,v(1)
704: PetscInt i,j,low,high,ii,ione,irow,row(1)
705: PetscScalar,pointer :: lx_v(:)
707: sc = solver%lambda
708: ione = 1
710: PetscCall(VecGetArrayReadF90(X1,lx_v,ierr))
711: PetscCall(VecGetOwnershipRange(X1,low,high,ierr))
713: ! Compute function over the locally owned part of the grid
714: ii = 0
715: do 20 irow=low,high-1
716: j = irow/solver%mx
717: i = mod(irow,solver%mx)
718: ii = ii + 1 ! one based local index
719: row(1) = irow
720: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
721: v(1) = 0.0
722: else
723: u = lx_v(ii)
724: v(1) = -sc*exp(u)
725: endif
726: PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr))
727: 20 continue
729: PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr))
731: PetscCall(VecAssemblyBegin(F1,ierr))
732: PetscCall(VecAssemblyEnd(F1,ierr))
734: ierr = 0
735: end subroutine FormFunctionNLTerm
737: !/*TEST
738: !
739: ! build:
740: ! requires: !single !complex
741: !
742: ! test:
743: ! nsize: 4
744: ! 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.
745: !
746: ! test:
747: ! args: -snes_linesearch_type {{l2 cp}separate output} -objective {{false true}shared output}
748: !
749: ! test:
750: ! args: -snes_linesearch_type bt -objective {{false true}separate output}
751: !
752: !TEST*/