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