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*/