Actual source code: ex1f.F90

petsc-master 2019-07-20
Report Typos and Errors
  1: !
  2: !
  3: !  Description: Uses the Newton method to solve a two-variable system.
  4: !
  5: !!/*T
  6: !  Concepts: SNES^basic uniprocessor example
  7: !  Processors: 1
  8: !T*/


 11: !
 12: ! -----------------------------------------------------------------------

 14:       program main
 15:  #include <petsc/finclude/petsc.h>
 16:       use petsc
 17:       implicit none


 20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 21: !                   Variable declarations
 22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 23: !
 24: !  Variables:
 25: !     snes        - nonlinear solver
 26: !     ksp        - linear solver
 27: !     pc          - preconditioner context
 28: !     ksp         - Krylov subspace method context
 29: !     x, r        - solution, residual vectors
 30: !     J           - Jacobian matrix
 31: !     its         - iterations for convergence
 32: !
 33:       SNES     snes
 34:       PC       pc
 35:       KSP      ksp
 36:       Vec      x,r
 37:       Mat      J
 38:       SNESLineSearch linesearch
 39:       PetscErrorCode  ierr
 40:       PetscInt its,i2,i20
 41:       PetscMPIInt size,rank
 42:       PetscScalar   pfive
 43:       PetscReal   tol
 44:       PetscBool   setls
 45:       PetscViewer viewer
 46:       double precision threshold,oldthreshold

 48: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 49: !  MUST be declared as external.

 51:       external FormFunction, FormJacobian, MyLineSearch

 53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54: !                   Macro definitions
 55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56: !
 57: !  Macros to make clearer the process of setting values in vectors and
 58: !  getting values from vectors.  These vectors are used in the routines
 59: !  FormFunction() and FormJacobian().
 60: !   - The element lx_a(ib) is element ib in the vector x
 61: !
 62: #define lx_a(ib) lx_v(lx_i + (ib))
 63: #define lf_a(ib) lf_v(lf_i + (ib))
 64: !
 65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66: !                 Beginning of program
 67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 69:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 70:       if (ierr .ne. 0) then
 71:         print*,'Unable to initialize PETSc'
 72:         stop
 73:       endif
 74:       call PetscLogNestedBegin(ierr);CHKERRA(ierr);
 75:       threshold = 1.0
 76:       call PetscLogSetThreshold(threshold,oldthreshold,ierr)
 77: ! dummy test of logging a reduction
 78:       PetscAReduce()
 79:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 80:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 81:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,1,'Uniprocessor example'); endif

 83:       i2  = 2
 84:       i20 = 20
 85: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 86: !  Create nonlinear solver context
 87: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 89:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92: !  Create matrix and vector data structures; set corresponding routines
 93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 95: !  Create vectors for solution and nonlinear function

 97:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
 98:       call VecDuplicate(x,r,ierr)

100: !  Create Jacobian matrix data structure

102:       call MatCreate(PETSC_COMM_SELF,J,ierr)
103:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
104:       call MatSetFromOptions(J,ierr)
105:       call MatSetUp(J,ierr)

107: !  Set function evaluation routine and vector

109:       call SNESSetFunction(snes,r,FormFunction,0,ierr)

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

113:       call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)

115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: !  Customize nonlinear solver; set runtime options
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

123:       call SNESGetKSP(snes,ksp,ierr)
124:       call KSPGetPC(ksp,pc,ierr)
125:       call PCSetType(pc,PCNONE,ierr)
126:       tol = 1.e-4
127:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                  &
128:      &                      PETSC_DEFAULT_REAL,i20,ierr)

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


137:       call SNESSetFromOptions(snes,ierr)

139:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
140:      &                         '-setls',setls,ierr)

142:       if (setls) then
143:         call SNESGetLineSearch(snes, linesearch, ierr)
144:         call SNESLineSearchSetType(linesearch, 'shell', ierr)
145:         call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,   &
146:      &                                      0, ierr)
147:       endif

149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !  Evaluate initial guess; then solve nonlinear system
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

158:       pfive = 0.5
159:       call VecSet(x,pfive,ierr)
160:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
161:       call SNESGetIterationNumber(snes,its,ierr);
162:       if (rank .eq. 0) then
163:          write(6,100) its
164:       endif
165:   100 format('Number of SNES iterations = ',i5)

167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: !  Free work space.  All PETSc objects should be destroyed when they
169: !  are no longer needed.
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

172:       call VecDestroy(x,ierr)
173:       call VecDestroy(r,ierr)
174:       call MatDestroy(J,ierr)
175:       call SNESDestroy(snes,ierr)
176:       call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)
177:       call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)
178:       call PetscLogView(viewer,ierr)
179:       call PetscViewerDestroy(viewer,ierr)
180:       call PetscFinalize(ierr)
181:       end
182: !
183: ! ------------------------------------------------------------------------
184: !
185: !  FormFunction - Evaluates nonlinear function, F(x).
186: !
187: !  Input Parameters:
188: !  snes - the SNES context
189: !  x - input vector
190: !  dummy - optional user-defined context (not used here)
191: !
192: !  Output Parameter:
193: !  f - function vector
194: !
195:       subroutine FormFunction(snes,x,f,dummy,ierr)
196:       use petscsnes
197:       implicit none

199:       SNES     snes
200:       Vec      x,f
201:       PetscErrorCode ierr
202:       integer dummy(*)

204: !  Declarations for use with local arrays

206:       PetscScalar  lx_v(2),lf_v(2)
207:       PetscOffset  lx_i,lf_i

209: !  Get pointers to vector data.
210: !    - For default PETSc vectors, VecGetArray() returns a pointer to
211: !      the data array.  Otherwise, the routine is implementation dependent.
212: !    - You MUST call VecRestoreArray() when you no longer need access to
213: !      the array.
214: !    - Note that the Fortran interface to VecGetArray() differs from the
215: !      C version.  See the Fortran chapter of the users manual for details.

217:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
218:       call VecGetArray(f,lf_v,lf_i,ierr)

220: !  Compute function

222:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
223:      &          + lx_a(1)*lx_a(2) - 3.0
224:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
225:      &          + lx_a(2)*lx_a(2) - 6.0

227: !  Restore vectors

229:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
230:       call VecRestoreArray(f,lf_v,lf_i,ierr)

232:       return
233:       end

235: ! ---------------------------------------------------------------------
236: !
237: !  FormJacobian - Evaluates Jacobian matrix.
238: !
239: !  Input Parameters:
240: !  snes - the SNES context
241: !  x - input vector
242: !  dummy - optional user-defined context (not used here)
243: !
244: !  Output Parameters:
245: !  A - Jacobian matrix
246: !  B - optionally different preconditioning matrix
247: !  flag - flag indicating matrix structure
248: !
249:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
250:       use petscsnes
251:       implicit none

253:       SNES         snes
254:       Vec          X
255:       Mat          jac,B
256:       PetscScalar  A(4)
257:       PetscErrorCode ierr
258:       PetscInt idx(2),i2
259:       integer dummy(*)

261: !  Declarations for use with local arrays

263:       PetscScalar lx_v(2)
264:       PetscOffset lx_i

266: !  Get pointer to vector data

268:       i2 = 2
269:       call VecGetArrayRead(x,lx_v,lx_i,ierr)

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

277:       idx(1) = 0
278:       idx(2) = 1
279:       A(1) = 2.0*lx_a(1) + lx_a(2)
280:       A(2) = lx_a(1)
281:       A(3) = lx_a(2)
282:       A(4) = lx_a(1) + 2.0*lx_a(2)
283:       call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)

285: !  Restore vector

287:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

289: !  Assemble matrix

291:       call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
292:       call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
293:       if (B .ne. jac) then
294:         call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
295:         call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
296:       endif

298:       return
299:       end


302:       subroutine MyLineSearch(linesearch, lctx, ierr)
303:       use petscsnes
304:       implicit none

306:       SNESLineSearch    linesearch
307:       SNES              snes
308:       integer           lctx
309:       Vec               x, f,g, y, w
310:       PetscReal         ynorm,gnorm,xnorm
311:       PetscBool         flag
312:       PetscErrorCode    ierr

314:       PetscScalar       mone

316:       mone = -1.0
317:       call SNESLineSearchGetSNES(linesearch, snes, ierr)
318:       call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
319:       call VecNorm(y,NORM_2,ynorm,ierr)
320:       call VecAXPY(x,mone,y,ierr)
321:       call SNESComputeFunction(snes,x,f,ierr)
322:       call VecNorm(f,NORM_2,gnorm,ierr)
323:       call VecNorm(x,NORM_2,xnorm,ierr)
324:       call VecNorm(y,NORM_2,ynorm,ierr)
325:       call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,      &
326:      & ierr)
327:       flag = PETSC_FALSE
328:       return
329:       end

331: !/*TEST
332: !
333: !   test:
334: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
335: !      requires: !single
336: !
337: !TEST*/