Actual source code: ex1f.F90
petsc-3.14.6 2021-03-30
1: !
2: !
3: ! Description: Uses the Newton method to solve a two-variable system.
4: !
5: !!/*T
6: ! Concepts: SNES^basic uniprocessor example
7: ! Processors: 1
8: !T*/
10: program main
11: #include <petsc/finclude/petsc.h>
12: use petsc
13: implicit none
15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16: ! Variable declarations
17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18: !
19: ! Variables:
20: ! snes - nonlinear solver
21: ! ksp - linear solver
22: ! pc - preconditioner context
23: ! ksp - Krylov subspace method context
24: ! x, r - solution, residual vectors
25: ! J - Jacobian matrix
26: ! its - iterations for convergence
27: !
28: SNES snes
29: PC pc
30: KSP ksp
31: Vec x,r
32: Mat J
33: SNESLineSearch linesearch
34: PetscErrorCode ierr
35: PetscInt its,i2,i20
36: PetscMPIInt size,rank
37: PetscScalar pfive
38: PetscReal tol
39: PetscBool setls
40: #if defined(PETSC_USE_LOG)
41: PetscViewer viewer
42: #endif
43: double precision threshold,oldthreshold
45: ! Note: Any user-defined Fortran routines (such as FormJacobian)
46: ! MUST be declared as external.
48: external FormFunction, FormJacobian, MyLineSearch
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Macro definitions
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: !
54: ! Macros to make clearer the process of setting values in vectors and
55: ! getting values from vectors. These vectors are used in the routines
56: ! FormFunction() and FormJacobian().
57: ! - The element lx_a(ib) is element ib in the vector x
58: !
59: #define lx_a(ib) lx_v(lx_i + (ib))
60: #define lf_a(ib) lf_v(lf_i + (ib))
61: !
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: ! Beginning of program
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
67: if (ierr .ne. 0) then
68: print*,'Unable to initialize PETSc'
69: stop
70: endif
71: call PetscLogNestedBegin(ierr);CHKERRA(ierr)
72: threshold = 1.0
73: call PetscLogSetThreshold(threshold,oldthreshold,ierr)
74: ! dummy test of logging a reduction
75: #if defined(PETSC_USE_LOG)
76: PetscAReduce()
77: #endif
78: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
79: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
80: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif
82: i2 = 2
83: i20 = 20
84: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
85: ! Create nonlinear solver context
86: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
88: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Create matrix and vector data structures; set corresponding routines
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: ! Create vectors for solution and nonlinear function
96: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
97: call VecDuplicate(x,r,ierr)
99: ! Create Jacobian matrix data structure
101: call MatCreate(PETSC_COMM_SELF,J,ierr)
102: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
103: call MatSetFromOptions(J,ierr)
104: call MatSetUp(J,ierr)
106: ! Set function evaluation routine and vector
108: call SNESSetFunction(snes,r,FormFunction,0,ierr)
110: ! Set Jacobian matrix data structure and Jacobian evaluation routine
112: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Customize nonlinear solver; set runtime options
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: ! Set linear solver defaults for this problem. By extracting the
119: ! KSP, KSP, and PC contexts from the SNES context, we can then
120: ! directly call any KSP, KSP, and PC routines to set various options.
122: call SNESGetKSP(snes,ksp,ierr)
123: call KSPGetPC(ksp,pc,ierr)
124: call PCSetType(pc,PCNONE,ierr)
125: tol = 1.e-4
126: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
127: & PETSC_DEFAULT_REAL,i20,ierr)
129: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
130: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
131: ! These options will override those specified above as long as
132: ! SNESSetFromOptions() is called _after_ any other customization
133: ! routines.
136: call SNESSetFromOptions(snes,ierr)
138: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
139: & '-setls',setls,ierr)
141: if (setls) then
142: call SNESGetLineSearch(snes, linesearch, ierr)
143: call SNESLineSearchSetType(linesearch, 'shell', ierr)
144: call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch, &
145: & 0, ierr)
146: endif
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: ! Evaluate initial guess; then solve nonlinear system
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: ! Note: The user should initialize the vector, x, with the initial guess
153: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
154: ! to employ an initial guess of zero, the user should explicitly set
155: ! this vector to zero by calling VecSet().
157: pfive = 0.5
158: call VecSet(x,pfive,ierr)
159: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
161: ! View solver converged reason; we could instead use the option -snes_converged_reason
162: call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr)
164: call SNESGetIterationNumber(snes,its,ierr);
165: if (rank .eq. 0) then
166: write(6,100) its
167: endif
168: 100 format('Number of SNES iterations = ',i5)
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: ! Free work space. All PETSc objects should be destroyed when they
172: ! are no longer needed.
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: call VecDestroy(x,ierr)
176: call VecDestroy(r,ierr)
177: call MatDestroy(J,ierr)
178: call SNESDestroy(snes,ierr)
179: #if defined(PETSC_USE_LOG)
180: call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)
181: call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)
182: call PetscLogView(viewer,ierr)
183: call PetscViewerDestroy(viewer,ierr)
184: #endif
185: call PetscFinalize(ierr)
186: end
187: !
188: ! ------------------------------------------------------------------------
189: !
190: ! FormFunction - Evaluates nonlinear function, F(x).
191: !
192: ! Input Parameters:
193: ! snes - the SNES context
194: ! x - input vector
195: ! dummy - optional user-defined context (not used here)
196: !
197: ! Output Parameter:
198: ! f - function vector
199: !
200: subroutine FormFunction(snes,x,f,dummy,ierr)
201: use petscsnes
202: implicit none
204: SNES snes
205: Vec x,f
206: PetscErrorCode ierr
207: integer dummy(*)
209: ! Declarations for use with local arrays
211: PetscScalar lx_v(2),lf_v(2)
212: PetscOffset lx_i,lf_i
214: ! Get pointers to vector data.
215: ! - For default PETSc vectors, VecGetArray() returns a pointer to
216: ! the data array. Otherwise, the routine is implementation dependent.
217: ! - You MUST call VecRestoreArray() when you no longer need access to
218: ! the array.
219: ! - Note that the Fortran interface to VecGetArray() differs from the
220: ! C version. See the Fortran chapter of the users manual for details.
222: call VecGetArrayRead(x,lx_v,lx_i,ierr)
223: call VecGetArray(f,lf_v,lf_i,ierr)
225: ! Compute function
227: lf_a(1) = lx_a(1)*lx_a(1) &
228: & + lx_a(1)*lx_a(2) - 3.0
229: lf_a(2) = lx_a(1)*lx_a(2) &
230: & + lx_a(2)*lx_a(2) - 6.0
232: ! Restore vectors
234: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
235: call VecRestoreArray(f,lf_v,lf_i,ierr)
237: return
238: end
240: ! ---------------------------------------------------------------------
241: !
242: ! FormJacobian - Evaluates Jacobian matrix.
243: !
244: ! Input Parameters:
245: ! snes - the SNES context
246: ! x - input vector
247: ! dummy - optional user-defined context (not used here)
248: !
249: ! Output Parameters:
250: ! A - Jacobian matrix
251: ! B - optionally different preconditioning matrix
252: ! flag - flag indicating matrix structure
253: !
254: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
255: use petscsnes
256: implicit none
258: SNES snes
259: Vec X
260: Mat jac,B
261: PetscScalar A(4)
262: PetscErrorCode ierr
263: PetscInt idx(2),i2
264: integer dummy(*)
266: ! Declarations for use with local arrays
268: PetscScalar lx_v(2)
269: PetscOffset lx_i
271: ! Get pointer to vector data
273: i2 = 2
274: call VecGetArrayRead(x,lx_v,lx_i,ierr)
276: ! Compute Jacobian entries and insert into matrix.
277: ! - Since this is such a small problem, we set all entries for
278: ! the matrix at once.
279: ! - Note that MatSetValues() uses 0-based row and column numbers
280: ! in Fortran as well as in C (as set here in the array idx).
282: idx(1) = 0
283: idx(2) = 1
284: A(1) = 2.0*lx_a(1) + lx_a(2)
285: A(2) = lx_a(1)
286: A(3) = lx_a(2)
287: A(4) = lx_a(1) + 2.0*lx_a(2)
288: call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
290: ! Restore vector
292: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
294: ! Assemble matrix
296: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
297: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
298: if (B .ne. jac) then
299: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
300: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
301: endif
303: return
304: end
307: subroutine MyLineSearch(linesearch, lctx, ierr)
308: use petscsnes
309: implicit none
311: SNESLineSearch linesearch
312: SNES snes
313: integer lctx
314: Vec x, f,g, y, w
315: PetscReal ynorm,gnorm,xnorm
316: PetscBool flag
317: PetscErrorCode ierr
319: PetscScalar mone
321: mone = -1.0
322: call SNESLineSearchGetSNES(linesearch, snes, ierr)
323: call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
324: call VecNorm(y,NORM_2,ynorm,ierr)
325: call VecAXPY(x,mone,y,ierr)
326: call SNESComputeFunction(snes,x,f,ierr)
327: call VecNorm(f,NORM_2,gnorm,ierr)
328: call VecNorm(x,NORM_2,xnorm,ierr)
329: call VecNorm(y,NORM_2,ynorm,ierr)
330: call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, &
331: & ierr)
332: flag = PETSC_FALSE
333: return
334: end
336: !/*TEST
337: !
338: ! test:
339: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
340: ! requires: !single
341: !
342: !TEST*/