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, intent(out) :: 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, intent(out) :: ierr
75: PetscInt idx(2)
76: integer dummy(*)
78: ! Declarations for use with local arrays
80: PetscScalar, pointer :: lx_v(:)
82: ! Get pointer to vector data
84: PetscCall(VecGetArrayRead(x, lx_v, ierr))
86: ! Compute Jacobian entries and insert into matrix.
87: ! - Since this is such a small problem, we set all entries for
88: ! the matrix at once.
89: ! - Note that MatSetValues() uses 0-based row and column numbers
90: ! in Fortran as well as in C (as set here in the array idx).
92: idx = [0, 1]
93: A = [2.0*lx_v(1) + lx_v(2), lx_v(1), lx_v(2), lx_v(1) + 2.0*lx_v(2)]
94: PetscCall(MatSetValues(B, 2_PETSC_INT_KIND, idx, 2_PETSC_INT_KIND, idx, A, INSERT_VALUES, ierr))
96: ! Restore vector
98: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
100: ! Assemble matrix
102: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
103: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
104: if (B /= jac) then
105: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
106: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
107: end if
109: end
111: subroutine MyLineSearch(linesearch, lctx, ierr)
113: SNESLineSearch linesearch
114: SNES snes
115: integer lctx
116: Vec x, f, g, y, w
117: PetscReal ynorm, gnorm, xnorm
118: PetscErrorCode, intent(out) :: ierr
120: PetscScalar, parameter :: mone = -1.0
122: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
123: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
124: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
125: PetscCall(VecAXPY(x, mone, y, ierr))
126: PetscCall(SNESComputeFunction(snes, x, f, ierr))
127: PetscCall(VecNorm(f, NORM_2, gnorm, ierr))
128: PetscCall(VecNorm(x, NORM_2, xnorm, ierr))
129: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
130: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, ierr))
131: end
132: end module
134: program main
135: use petsc
136: use ex1fmodule
137: implicit none
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: ! Variable declarations
141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: !
143: ! Variables:
144: ! snes - nonlinear solver
145: ! ksp - linear solver
146: ! pc - preconditioner context
147: ! ksp - Krylov subspace method context
148: ! x, r - solution, residual vectors
149: ! J - Jacobian matrix
150: ! its - iterations for convergence
151: !
152: SNES snes
153: PC pc
154: KSP ksp
155: Vec x, r
156: Mat J
157: SNESLineSearch linesearch
158: PetscErrorCode ierr
159: PetscInt its
160: PetscMPIInt size, rank
161: PetscScalar, parameter :: pfive = 0.5
162: PetscReal, parameter :: tol = 1.e-4
163: PetscBool setls
164: PetscReal, pointer :: rhistory(:)
165: PetscInt, pointer :: itshistory(:)
166: PetscInt nhistory
167: #if defined(PETSC_USE_LOG)
168: PetscViewer viewer
169: #endif
170: double precision threshold, oldthreshold
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: ! Beginning of program
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: PetscCallA(PetscInitialize(ierr))
177: PetscCallA(PetscLogNestedBegin(ierr))
178: threshold = 1.0
179: PetscCallA(PetscLogSetThreshold(threshold, oldthreshold, ierr))
180: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
181: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
182: PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'Uniprocessor example')
184: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! Create nonlinear solver context
186: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
188: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
190: PetscCallA(SNESSetConvergenceHistory(snes, PETSC_NULL_REAL_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_DECIDE, PETSC_FALSE, ierr))
192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: ! Create matrix and vector data structures; set corresponding routines
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: ! Create vectors for solution and nonlinear function
198: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, 2_PETSC_INT_KIND, x, ierr))
199: PetscCallA(VecDuplicate(x, r, ierr))
201: ! Create Jacobian matrix data structure
203: PetscCallA(MatCreate(PETSC_COMM_SELF, J, ierr))
204: PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
205: PetscCallA(MatSetFromOptions(J, ierr))
206: PetscCallA(MatSetUp(J, ierr))
208: ! Set function evaluation routine and vector
210: PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))
212: ! Set Jacobian matrix data structure and Jacobian evaluation routine
214: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: ! Customize nonlinear solver; set runtime options
218: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: ! Set linear solver defaults for this problem. By extracting the
221: ! KSP, KSP, and PC contexts from the SNES context, we can then
222: ! directly call any KSP, KSP, and PC routines to set various options.
224: PetscCallA(SNESGetKSP(snes, ksp, ierr))
225: PetscCallA(KSPGetPC(ksp, pc, ierr))
226: PetscCallA(PCSetType(pc, PCNONE, ierr))
227: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, 20_PETSC_INT_KIND, ierr))
229: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
230: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
231: ! These options will override those specified above as long as
232: ! SNESSetFromOptions() is called _after_ any other customization
233: ! routines.
235: PetscCallA(SNESSetFromOptions(snes, ierr))
237: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-setls', setls, ierr))
239: if (setls) then
240: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
241: PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
242: PetscCallA(SNESLineSearchShellSetApply(linesearch, MyLineSearch, 0, ierr))
243: end if
245: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
246: ! Evaluate initial guess; then solve nonlinear system
247: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249: ! Note: The user should initialize the vector, x, with the initial guess
250: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
251: ! to employ an initial guess of zero, the user should explicitly set
252: ! this vector to zero by calling VecSet().
254: PetscCallA(VecSet(x, pfive, ierr))
255: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
257: PetscCallA(SNESGetConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
258: PetscCallA(SNESRestoreConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
260: ! View solver converged reason; we could instead use the option -snes_converged_reason
261: PetscCallA(SNESConvergedReasonView(snes, PETSC_VIEWER_STDOUT_WORLD, ierr))
263: PetscCallA(SNESGetIterationNumber(snes, its, ierr))
264: if (rank == 0) then
265: write (6, 100) its
266: end if
267: 100 format('Number of SNES iterations = ', i5)
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: ! Free work space. All PETSc objects should be destroyed when they
271: ! are no longer needed.
272: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: PetscCallA(VecDestroy(x, ierr))
275: PetscCallA(VecDestroy(r, ierr))
276: PetscCallA(MatDestroy(J, ierr))
277: PetscCallA(SNESDestroy(snes, ierr))
278: #if defined(PETSC_USE_LOG)
279: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD, 'filename.xml', viewer, ierr))
280: PetscCallA(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_XML, ierr))
281: PetscCallA(PetscLogView(viewer, ierr))
282: PetscCallA(PetscViewerDestroy(viewer, ierr))
283: #endif
284: PetscCallA(PetscFinalize(ierr))
285: end
286: !/*TEST
287: !
288: ! test:
289: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
290: ! requires: !single
291: !
292: !TEST*/