Actual source code: ex1f.F90
1: !
2: ! Description: Uses the Newton method to solve a two-variable system.
3: !
4: #include <petsc/finclude/petsc.h>
5: module ex1fmodule
6: use petscvec
7: use petscsnesdef
8: use petscvec
9: use petscmat
10: implicit none
12: contains
13: !
14: ! ------------------------------------------------------------------------
15: !
16: ! FormFunction - Evaluates nonlinear function, F(x).
17: !
18: ! Input Parameters:
19: ! snes - the SNES context
20: ! x - input vector
21: ! dummy - optional user-defined context (not used here)
22: !
23: ! Output Parameter:
24: ! f - function vector
25: !
26: subroutine FormFunction(snes, x, f, dummy, ierr)
27: SNES snes
28: Vec x, f
29: PetscErrorCode ierr
30: integer dummy(*)
32: ! Declarations for use with local arrays
33: PetscScalar, pointer :: lx_v(:), lf_v(:)
35: ! Get pointers to vector data.
36: ! - VecGetArray() returns a pointer to the data array.
37: ! - You MUST call VecRestoreArray() when you no longer need access to
38: ! the array.
40: PetscCall(VecGetArrayRead(x, lx_v, ierr))
41: PetscCall(VecGetArray(f, lf_v, ierr))
43: ! Compute function
45: lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
46: lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0
48: ! Restore vectors
50: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
51: PetscCall(VecRestoreArray(f, lf_v, ierr))
53: end
55: ! ---------------------------------------------------------------------
56: !
57: ! FormJacobian - Evaluates Jacobian matrix.
58: !
59: ! Input Parameters:
60: ! snes - the SNES context
61: ! x - input vector
62: ! dummy - optional user-defined context (not used here)
63: !
64: ! Output Parameters:
65: ! A - Jacobian matrix
66: ! B - optionally different matrix used to construct the preconditioner
67: !
68: subroutine FormJacobian(snes, X, jac, B, dummy, ierr)
70: SNES snes
71: Vec X
72: Mat jac, B
73: PetscScalar A(4)
74: PetscErrorCode ierr
75: PetscInt idx(2), i2
76: integer dummy(*)
78: ! Declarations for use with local arrays
80: PetscScalar, pointer :: lx_v(:)
82: ! Get pointer to vector data
84: i2 = 2
85: PetscCall(VecGetArrayRead(x, lx_v, ierr))
87: ! Compute Jacobian entries and insert into matrix.
88: ! - Since this is such a small problem, we set all entries for
89: ! the matrix at once.
90: ! - Note that MatSetValues() uses 0-based row and column numbers
91: ! in Fortran as well as in C (as set here in the array idx).
93: idx(1) = 0
94: idx(2) = 1
95: A(1) = 2.0*lx_v(1) + lx_v(2)
96: A(2) = lx_v(1)
97: A(3) = lx_v(2)
98: A(4) = lx_v(1) + 2.0*lx_v(2)
99: PetscCall(MatSetValues(B, i2, idx, i2, idx, A, INSERT_VALUES, ierr))
101: ! Restore vector
103: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
105: ! Assemble matrix
107: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
108: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
109: if (B /= jac) then
110: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
111: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
112: end if
114: end
116: subroutine MyLineSearch(linesearch, lctx, ierr)
118: SNESLineSearch linesearch
119: SNES snes
120: integer lctx
121: Vec x, f, g, y, w
122: PetscReal ynorm, gnorm, xnorm
123: PetscErrorCode ierr
125: PetscScalar mone
127: mone = -1.0
128: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
129: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
130: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
131: PetscCall(VecAXPY(x, mone, y, ierr))
132: PetscCall(SNESComputeFunction(snes, x, f, ierr))
133: PetscCall(VecNorm(f, NORM_2, gnorm, ierr))
134: PetscCall(VecNorm(x, NORM_2, xnorm, ierr))
135: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
136: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, ierr))
137: end
138: end module
140: program main
141: use petsc
142: use ex1fmodule
143: implicit none
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: ! Variable declarations
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: !
149: ! Variables:
150: ! snes - nonlinear solver
151: ! ksp - linear solver
152: ! pc - preconditioner context
153: ! ksp - Krylov subspace method context
154: ! x, r - solution, residual vectors
155: ! J - Jacobian matrix
156: ! its - iterations for convergence
157: !
158: SNES snes
159: PC pc
160: KSP ksp
161: Vec x, r
162: Mat J
163: SNESLineSearch linesearch
164: PetscErrorCode ierr
165: PetscInt its, i2, i20
166: PetscMPIInt size, rank
167: PetscScalar pfive
168: PetscReal tol
169: PetscBool setls
170: PetscReal, pointer :: rhistory(:)
171: PetscInt, pointer :: itshistory(:)
172: PetscInt nhistory
173: #if defined(PETSC_USE_LOG)
174: PetscViewer viewer
175: #endif
176: double precision threshold, oldthreshold
178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: ! Beginning of program
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: PetscCallA(PetscInitialize(ierr))
183: PetscCallA(PetscLogNestedBegin(ierr))
184: threshold = 1.0
185: PetscCallA(PetscLogSetThreshold(threshold, oldthreshold, ierr))
186: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
187: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
188: PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'Uniprocessor example')
190: i2 = 2
191: i20 = 20
192: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
193: ! Create nonlinear solver context
194: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
196: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
198: PetscCallA(SNESSetConvergenceHistory(snes, PETSC_NULL_REAL_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_DECIDE, PETSC_FALSE, ierr))
200: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: ! Create matrix and vector data structures; set corresponding routines
202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: ! Create vectors for solution and nonlinear function
206: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i2, x, ierr))
207: PetscCallA(VecDuplicate(x, r, ierr))
209: ! Create Jacobian matrix data structure
211: PetscCallA(MatCreate(PETSC_COMM_SELF, J, ierr))
212: PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, i2, i2, ierr))
213: PetscCallA(MatSetFromOptions(J, ierr))
214: PetscCallA(MatSetUp(J, ierr))
216: ! Set function evaluation routine and vector
218: PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))
220: ! Set Jacobian matrix data structure and Jacobian evaluation routine
222: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
224: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225: ! Customize nonlinear solver; set runtime options
226: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: ! Set linear solver defaults for this problem. By extracting the
229: ! KSP, KSP, and PC contexts from the SNES context, we can then
230: ! directly call any KSP, KSP, and PC routines to set various options.
232: PetscCallA(SNESGetKSP(snes, ksp, ierr))
233: PetscCallA(KSPGetPC(ksp, pc, ierr))
234: PetscCallA(PCSetType(pc, PCNONE, ierr))
235: tol = 1.e-4
236: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, i20, ierr))
238: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
239: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
240: ! These options will override those specified above as long as
241: ! SNESSetFromOptions() is called _after_ any other customization
242: ! routines.
244: PetscCallA(SNESSetFromOptions(snes, ierr))
246: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-setls', setls, ierr))
248: if (setls) then
249: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
250: PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
251: PetscCallA(SNESLineSearchShellSetApply(linesearch, MyLineSearch, 0, ierr))
252: end if
254: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: ! Evaluate initial guess; then solve nonlinear system
256: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: ! Note: The user should initialize the vector, x, with the initial guess
259: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
260: ! to employ an initial guess of zero, the user should explicitly set
261: ! this vector to zero by calling VecSet().
263: pfive = 0.5
264: PetscCallA(VecSet(x, pfive, ierr))
265: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
267: PetscCallA(SNESGetConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
268: PetscCallA(SNESRestoreConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
270: ! View solver converged reason; we could instead use the option -snes_converged_reason
271: PetscCallA(SNESConvergedReasonView(snes, PETSC_VIEWER_STDOUT_WORLD, ierr))
273: PetscCallA(SNESGetIterationNumber(snes, its, ierr))
274: if (rank == 0) then
275: write (6, 100) its
276: end if
277: 100 format('Number of SNES iterations = ', i5)
279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: ! Free work space. All PETSc objects should be destroyed when they
281: ! are no longer needed.
282: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284: PetscCallA(VecDestroy(x, ierr))
285: PetscCallA(VecDestroy(r, ierr))
286: PetscCallA(MatDestroy(J, ierr))
287: PetscCallA(SNESDestroy(snes, ierr))
288: #if defined(PETSC_USE_LOG)
289: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD, 'filename.xml', viewer, ierr))
290: PetscCallA(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_XML, ierr))
291: PetscCallA(PetscLogView(viewer, ierr))
292: PetscCallA(PetscViewerDestroy(viewer, ierr))
293: #endif
294: PetscCallA(PetscFinalize(ierr))
295: end
296: !/*TEST
297: !
298: ! test:
299: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
300: ! requires: !single
301: !
302: !TEST*/