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: ! dummy test of logging a reduction
 55: #if defined(PETSC_USE_LOG)
 56:       ierr = PetscAReduce()
 57: #endif
 58:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 59:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 60:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif

 62:       i2  = 2
 63:       i20 = 20
 64: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 65: !  Create nonlinear solver context
 66: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 68:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

 70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71: !  Create matrix and vector data structures; set corresponding routines
 72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 74: !  Create vectors for solution and nonlinear function

 76:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr))
 77:       PetscCallA(VecDuplicate(x,r,ierr))

 79: !  Create Jacobian matrix data structure

 81:       PetscCallA(MatCreate(PETSC_COMM_SELF,J,ierr))
 82:       PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr))
 83:       PetscCallA(MatSetFromOptions(J,ierr))
 84:       PetscCallA(MatSetUp(J,ierr))

 86: !  Set function evaluation routine and vector

 88:       PetscCallA(SNESSetFunction(snes,r,FormFunction,0,ierr))

 90: !  Set Jacobian matrix data structure and Jacobian evaluation routine

 92:       PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))

 94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95: !  Customize nonlinear solver; set runtime options
 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 98: !  Set linear solver defaults for this problem. By extracting the
 99: !  KSP, KSP, and PC contexts from the SNES context, we can then
100: !  directly call any KSP, KSP, and PC routines to set various options.

102:       PetscCallA(SNESGetKSP(snes,ksp,ierr))
103:       PetscCallA(KSPGetPC(ksp,pc,ierr))
104:       PetscCallA(PCSetType(pc,PCNONE,ierr))
105:       tol = 1.e-4
106:       PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,i20,ierr))

108: !  Set SNES/KSP/KSP/PC runtime options, e.g.,
109: !      -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
110: !  These options will override those specified above as long as
111: !  SNESSetFromOptions() is called _after_ any other customization
112: !  routines.

114:       PetscCallA(SNESSetFromOptions(snes,ierr))

116:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-setls',setls,ierr))

118:       if (setls) then
119:         PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
120:         PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
121:         PetscCallA(SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,0,ierr))
122:       endif

124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: !  Evaluate initial guess; then solve nonlinear system
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

128: !  Note: The user should initialize the vector, x, with the initial guess
129: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
130: !  to employ an initial guess of zero, the user should explicitly set
131: !  this vector to zero by calling VecSet().

133:       pfive = 0.5
134:       PetscCallA(VecSet(x,pfive,ierr))
135:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))

137: !  View solver converged reason; we could instead use the option -snes_converged_reason
138:       PetscCallA(SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr))

140:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))
141:       if (rank .eq. 0) then
142:          write(6,100) its
143:       endif
144:   100 format('Number of SNES iterations = ',i5)

146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: !  Free work space.  All PETSc objects should be destroyed when they
148: !  are no longer needed.
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

151:       PetscCallA(VecDestroy(x,ierr))
152:       PetscCallA(VecDestroy(r,ierr))
153:       PetscCallA(MatDestroy(J,ierr))
154:       PetscCallA(SNESDestroy(snes,ierr))
155: #if defined(PETSC_USE_LOG)
156:       PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr))
157:       PetscCallA(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr))
158:       PetscCallA(PetscLogView(viewer,ierr))
159:       PetscCallA(PetscViewerDestroy(viewer,ierr))
160: #endif
161:       PetscCallA(PetscFinalize(ierr))
162:       end
163: !
164: ! ------------------------------------------------------------------------
165: !
166: !  FormFunction - Evaluates nonlinear function, F(x).
167: !
168: !  Input Parameters:
169: !  snes - the SNES context
170: !  x - input vector
171: !  dummy - optional user-defined context (not used here)
172: !
173: !  Output Parameter:
174: !  f - function vector
175: !
176:       subroutine FormFunction(snes,x,f,dummy,ierr)
177:       use petscsnes
178:       implicit none

180:       SNES     snes
181:       Vec      x,f
182:       PetscErrorCode ierr
183:       integer dummy(*)

185: !  Declarations for use with local arrays
186:       PetscScalar,pointer :: lx_v(:),lf_v(:)

188: !  Get pointers to vector data.
189: !    - VecGetArrayF90() returns a pointer to the data array.
190: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
191: !      the array.

193:       PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
194:       PetscCall(VecGetArrayF90(f,lf_v,ierr))

196: !  Compute function

198:       lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
199:       lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0

201: !  Restore vectors

203:       PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
204:       PetscCall(VecRestoreArrayF90(f,lf_v,ierr))

206:       return
207:       end

209: ! ---------------------------------------------------------------------
210: !
211: !  FormJacobian - Evaluates Jacobian matrix.
212: !
213: !  Input Parameters:
214: !  snes - the SNES context
215: !  x - input vector
216: !  dummy - optional user-defined context (not used here)
217: !
218: !  Output Parameters:
219: !  A - Jacobian matrix
220: !  B - optionally different preconditioning matrix
221: !
222:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
223:       use petscsnes
224:       implicit none

226:       SNES         snes
227:       Vec          X
228:       Mat          jac,B
229:       PetscScalar  A(4)
230:       PetscErrorCode ierr
231:       PetscInt idx(2),i2
232:       integer dummy(*)

234: !  Declarations for use with local arrays

236:       PetscScalar,pointer :: lx_v(:)

238: !  Get pointer to vector data

240:       i2 = 2
241:       PetscCall(VecGetArrayReadF90(x,lx_v,ierr))

243: !  Compute Jacobian entries and insert into matrix.
244: !   - Since this is such a small problem, we set all entries for
245: !     the matrix at once.
246: !   - Note that MatSetValues() uses 0-based row and column numbers
247: !     in Fortran as well as in C (as set here in the array idx).

249:       idx(1) = 0
250:       idx(2) = 1
251:       A(1) = 2.0*lx_v(1) + lx_v(2)
252:       A(2) = lx_v(1)
253:       A(3) = lx_v(2)
254:       A(4) = lx_v(1) + 2.0*lx_v(2)
255:       PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr))

257: !  Restore vector

259:       PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))

261: !  Assemble matrix

263:       PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
264:       PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
265:       if (B .ne. jac) then
266:         PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
267:         PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
268:       endif

270:       return
271:       end

273:       subroutine MyLineSearch(linesearch, lctx, ierr)
274:       use petscsnes
275:       implicit none

277:       SNESLineSearch    linesearch
278:       SNES              snes
279:       integer           lctx
280:       Vec               x, f,g, y, w
281:       PetscReal         ynorm,gnorm,xnorm
282:       PetscErrorCode    ierr

284:       PetscScalar       mone

286:       mone = -1.0
287:       PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
288:       PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
289:       PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
290:       PetscCall(VecAXPY(x,mone,y,ierr))
291:       PetscCall(SNESComputeFunction(snes,x,f,ierr))
292:       PetscCall(VecNorm(f,NORM_2,gnorm,ierr))
293:       PetscCall(VecNorm(x,NORM_2,xnorm,ierr))
294:       PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
295:       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr))
296:       return
297:       end

299: !/*TEST
300: !
301: !   test:
302: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
303: !      requires: !single
304: !
305: !TEST*/