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