Actual source code: ex1f.F90
petsc-3.8.4 2018-03-24
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*/
9: !
10: ! -----------------------------------------------------------------------
12: program main
13: #include <petsc/finclude/petsc.h>
14: use petsc
15: implicit none
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: ! Variable declarations
20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
21: !
22: ! Variables:
23: ! snes - nonlinear solver
24: ! ksp - linear solver
25: ! pc - preconditioner context
26: ! ksp - Krylov subspace method context
27: ! x, r - solution, residual vectors
28: ! J - Jacobian matrix
29: ! its - iterations for convergence
30: !
31: SNES snes
32: PC pc
33: KSP ksp
34: Vec x,r
35: Mat J
36: SNESLineSearch linesearch
37: PetscErrorCode ierr
38: PetscInt its,i2,i20
39: PetscMPIInt size,rank
40: PetscScalar pfive
41: PetscReal tol
42: PetscBool setls
44: ! Note: Any user-defined Fortran routines (such as FormJacobian)
45: ! MUST be declared as external.
47: external FormFunction, FormJacobian, MyLineSearch
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Macro definitions
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: !
53: ! Macros to make clearer the process of setting values in vectors and
54: ! getting values from vectors. These vectors are used in the routines
55: ! FormFunction() and FormJacobian().
56: ! - The element lx_a(ib) is element ib in the vector x
57: !
58: #define lx_a(ib) lx_v(lx_i + (ib))
59: #define lf_a(ib) lf_v(lf_i + (ib))
60: !
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: ! Beginning of program
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
66: if (ierr .ne. 0) then
67: print*,'Unable to initialize PETSc'
68: stop
69: endif
70: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
71: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
72: if (size .ne. 1) then
73: if (rank .eq. 0) then
74: write(6,*) 'This is a uniprocessor example only!'
75: endif
76: SETERRA(PETSC_COMM_SELF,1,'')
77: endif
79: i2 = 2
80: i20 = 20
81: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
82: ! Create nonlinear solver context
83: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
85: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: ! Create matrix and vector data structures; set corresponding routines
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Create vectors for solution and nonlinear function
93: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
94: call VecDuplicate(x,r,ierr)
96: ! Create Jacobian matrix data structure
98: call MatCreate(PETSC_COMM_SELF,J,ierr)
99: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
100: call MatSetFromOptions(J,ierr)
101: call MatSetUp(J,ierr)
103: ! Set function evaluation routine and vector
105: call SNESSetFunction(snes,r,FormFunction,0,ierr)
107: ! Set Jacobian matrix data structure and Jacobian evaluation routine
109: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: ! Customize nonlinear solver; set runtime options
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Set linear solver defaults for this problem. By extracting the
116: ! KSP, KSP, and PC contexts from the SNES context, we can then
117: ! directly call any KSP, KSP, and PC routines to set various options.
119: call SNESGetKSP(snes,ksp,ierr)
120: call KSPGetPC(ksp,pc,ierr)
121: call PCSetType(pc,PCNONE,ierr)
122: tol = 1.e-4
123: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
124: & PETSC_DEFAULT_REAL,i20,ierr)
126: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
127: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
128: ! These options will override those specified above as long as
129: ! SNESSetFromOptions() is called _after_ any other customization
130: ! routines.
133: call SNESSetFromOptions(snes,ierr)
135: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
136: & '-setls',setls,ierr)
138: if (setls) then
139: call SNESGetLineSearch(snes, linesearch, ierr)
140: call SNESLineSearchSetType(linesearch, 'shell', ierr)
141: call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch, &
142: & 0, ierr)
143: endif
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: ! Evaluate initial guess; then solve nonlinear system
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: ! Note: The user should initialize the vector, x, with the initial guess
150: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
151: ! to employ an initial guess of zero, the user should explicitly set
152: ! this vector to zero by calling VecSet().
154: pfive = 0.5
155: call VecSet(x,pfive,ierr)
156: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
157: call SNESGetIterationNumber(snes,its,ierr);
158: if (rank .eq. 0) then
159: write(6,100) its
160: endif
161: 100 format('Number of SNES iterations = ',i5)
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Free work space. All PETSc objects should be destroyed when they
165: ! are no longer needed.
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: call VecDestroy(x,ierr)
169: call VecDestroy(r,ierr)
170: call MatDestroy(J,ierr)
171: call SNESDestroy(snes,ierr)
172: call PetscFinalize(ierr)
173: end
174: !
175: ! ------------------------------------------------------------------------
176: !
177: ! FormFunction - Evaluates nonlinear function, F(x).
178: !
179: ! Input Parameters:
180: ! snes - the SNES context
181: ! x - input vector
182: ! dummy - optional user-defined context (not used here)
183: !
184: ! Output Parameter:
185: ! f - function vector
186: !
187: subroutine FormFunction(snes,x,f,dummy,ierr)
188: use petscsnes
189: implicit none
191: SNES snes
192: Vec x,f
193: PetscErrorCode ierr
194: integer dummy(*)
196: ! Declarations for use with local arrays
198: PetscScalar lx_v(2),lf_v(2)
199: PetscOffset lx_i,lf_i
201: ! Get pointers to vector data.
202: ! - For default PETSc vectors, VecGetArray() returns a pointer to
203: ! the data array. Otherwise, the routine is implementation dependent.
204: ! - You MUST call VecRestoreArray() when you no longer need access to
205: ! the array.
206: ! - Note that the Fortran interface to VecGetArray() differs from the
207: ! C version. See the Fortran chapter of the users manual for details.
209: call VecGetArrayRead(x,lx_v,lx_i,ierr)
210: call VecGetArray(f,lf_v,lf_i,ierr)
212: ! Compute function
214: lf_a(1) = lx_a(1)*lx_a(1) &
215: & + lx_a(1)*lx_a(2) - 3.0
216: lf_a(2) = lx_a(1)*lx_a(2) &
217: & + lx_a(2)*lx_a(2) - 6.0
219: ! Restore vectors
221: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
222: call VecRestoreArray(f,lf_v,lf_i,ierr)
224: return
225: end
227: ! ---------------------------------------------------------------------
228: !
229: ! FormJacobian - Evaluates Jacobian matrix.
230: !
231: ! Input Parameters:
232: ! snes - the SNES context
233: ! x - input vector
234: ! dummy - optional user-defined context (not used here)
235: !
236: ! Output Parameters:
237: ! A - Jacobian matrix
238: ! B - optionally different preconditioning matrix
239: ! flag - flag indicating matrix structure
240: !
241: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
242: use petscsnes
243: implicit none
245: SNES snes
246: Vec X
247: Mat jac,B
248: PetscScalar A(4)
249: PetscErrorCode ierr
250: PetscInt idx(2),i2
251: integer dummy(*)
253: ! Declarations for use with local arrays
255: PetscScalar lx_v(2)
256: PetscOffset lx_i
258: ! Get pointer to vector data
260: i2 = 2
261: call VecGetArrayRead(x,lx_v,lx_i,ierr)
263: ! Compute Jacobian entries and insert into matrix.
264: ! - Since this is such a small problem, we set all entries for
265: ! the matrix at once.
266: ! - Note that MatSetValues() uses 0-based row and column numbers
267: ! in Fortran as well as in C (as set here in the array idx).
269: idx(1) = 0
270: idx(2) = 1
271: A(1) = 2.0*lx_a(1) + lx_a(2)
272: A(2) = lx_a(1)
273: A(3) = lx_a(2)
274: A(4) = lx_a(1) + 2.0*lx_a(2)
275: call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
277: ! Restore vector
279: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
281: ! Assemble matrix
283: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
284: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
285: if (B .ne. jac) then
286: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
287: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
288: endif
290: return
291: end
294: subroutine MyLineSearch(linesearch, lctx, ierr)
295: use petscsnes
296: implicit none
298: SNESLineSearch linesearch
299: SNES snes
300: integer lctx
301: Vec x, f,g, y, w
302: PetscReal ynorm,gnorm,xnorm
303: PetscBool flag
304: PetscErrorCode ierr
306: PetscScalar mone
308: mone = -1.0
309: call SNESLineSearchGetSNES(linesearch, snes, ierr)
310: call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
311: call VecNorm(y,NORM_2,ynorm,ierr)
312: call VecAXPY(x,mone,y,ierr)
313: call SNESComputeFunction(snes,x,f,ierr)
314: call VecNorm(f,NORM_2,gnorm,ierr)
315: call VecNorm(x,NORM_2,xnorm,ierr)
316: call VecNorm(y,NORM_2,ynorm,ierr)
317: call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, &
318: & ierr)
319: flag = PETSC_FALSE
320: return
321: end