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:       return
203:       end

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

222:       SNES         snes
223:       Vec          X
224:       Mat          jac,B
225:       PetscScalar  A(4)
226:       PetscErrorCode ierr
227:       PetscInt idx(2),i2
228:       integer dummy(*)

230: !  Declarations for use with local arrays

232:       PetscScalar,pointer :: lx_v(:)

234: !  Get pointer to vector data

236:       i2 = 2
237:       PetscCall(VecGetArrayReadF90(x,lx_v,ierr))

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

245:       idx(1) = 0
246:       idx(2) = 1
247:       A(1) = 2.0*lx_v(1) + lx_v(2)
248:       A(2) = lx_v(1)
249:       A(3) = lx_v(2)
250:       A(4) = lx_v(1) + 2.0*lx_v(2)
251:       PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr))

253: !  Restore vector

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

257: !  Assemble matrix

259:       PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
260:       PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
261:       if (B .ne. jac) then
262:         PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
263:         PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
264:       endif

266:       return
267:       end

269:       subroutine MyLineSearch(linesearch, lctx, ierr)
270:       use petscsnes
271:       implicit none

273:       SNESLineSearch    linesearch
274:       SNES              snes
275:       integer           lctx
276:       Vec               x, f,g, y, w
277:       PetscReal         ynorm,gnorm,xnorm
278:       PetscErrorCode    ierr

280:       PetscScalar       mone

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

295: !/*TEST
296: !
297: !   test:
298: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
299: !      requires: !single
300: !
301: !TEST*/