Actual source code: eptorsion2f.F90

  1: !  Program usage: mpiexec -n <proc> eptorsion2f [all TAO options]
  2: !
  3: !  Description:  This example demonstrates use of the TAO package to solve
  4: !  unconstrained minimization problems in parallel.  This example is based
  5: !  on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
  6: !  The command line options are:
  7: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
  8: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
  9: !    -par <param>, where <param> = angle of twist per unit length
 10: !
 11: !
 12: ! ----------------------------------------------------------------------
 13: !
 14: !  Elastic-plastic torsion problem.
 15: !
 16: !  The elastic plastic torsion problem arises from the deconverged
 17: !  of the stress field on an infinitely long cylindrical bar, which is
 18: !  equivalent to the solution of the following problem:
 19: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
 20: !  where C is the torsion angle per unit length.
 21: !
 22: !  The C version of this code is eptorsion2.c
 23: !
 24: ! ----------------------------------------------------------------------

 26:       module eptorsion2fmodule
 27: #include "petsc/finclude/petsctao.h"
 28:       use petscdmda
 29:       use petsctao
 30:       implicit none

 32:       Vec              localX
 33:       DM               dm
 34:       PetscReal      param
 35:       PetscInt         mx, my
 36:       end module

 38:       use eptorsion2fmodule
 39:       implicit none
 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !                   Variable declarations
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !
 44: !  See additional variable declarations in the file eptorsion2f.h
 45: !
 46:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
 47:       Vec              x              ! solution vector
 48:       Mat              H              ! hessian matrix
 49:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 50:       Tao        tao            ! Tao solver context
 51:       PetscBool        flg
 52:       PetscInt         i1
 53:       PetscInt         dummy

 55: !  Note: Any user-defined Fortran routines (such as FormGradient)
 56: !  MUST be declared as external.

 58:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
 59:       external Monitor,ConvergenceTest

 61:       i1 = 1

 63: !     Initialize TAO, PETSc  contexts
 64:       PetscCallA(PetscInitialize(ierr))

 66: !     Specify default parameters
 67:       param = 5.0
 68:       mx = 10
 69:       my = 10
 70:       Nx = PETSC_DECIDE
 71:       Ny = PETSC_DECIDE

 73: !     Check for any command line arguments that might override defaults
 74:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
 75:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
 76:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',param,flg,ierr))

 78: !     Set up distributed array and vectors
 79:       PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,dm,ierr))
 80:       PetscCallA(DMSetFromOptions(dm,ierr))
 81:       PetscCallA(DMSetUp(dm,ierr))

 83: !     Create vectors
 84:       PetscCallA(DMCreateGlobalVector(dm,x,ierr))
 85:       PetscCallA(DMCreateLocalVector(dm,localX,ierr))

 87: !     Create Hessian
 88:       PetscCallA(DMCreateMatrix(dm,H,ierr))
 89:       PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

 91: !     The TAO code begins here

 93: !     Create TAO solver
 94:       PetscCallA(TaoCreate(PETSC_COMM_WORLD,tao,ierr))
 95:       PetscCallA(TaoSetType(tao,TAOCG,ierr))

 97: !     Set routines for function and gradient evaluation

 99:       PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
100:       PetscCallA(TaoSetHessian(tao,H,H,ComputeHessian,0,ierr))

102: !     Set initial guess
103:       PetscCallA(FormInitialGuess(x,ierr))
104:       PetscCallA(TaoSetSolution(tao,x,ierr))

106:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testmonitor',flg,ierr))
107:       if (flg) then
108:          PetscCallA(TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,ierr))
109:       endif

111:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testconvergence',flg, ierr))
112:       if (flg) then
113:          PetscCallA(TaoSetConvergenceTest(tao,ConvergenceTest,dummy,ierr))
114:       endif

116: !     Check for any TAO command line options
117:       PetscCallA(TaoSetFromOptions(tao,ierr))

119: !     SOLVE THE APPLICATION
120:       PetscCallA(TaoSolve(tao,ierr))

122: !     Free TAO data structures
123:       PetscCallA(TaoDestroy(tao,ierr))

125: !     Free PETSc data structures
126:       PetscCallA(VecDestroy(x,ierr))
127:       PetscCallA(VecDestroy(localX,ierr))
128:       PetscCallA(MatDestroy(H,ierr))
129:       PetscCallA(DMDestroy(dm,ierr))

131: !     Finalize TAO and PETSc
132:       PetscCallA(PetscFinalize(ierr))
133:       end

135: ! ---------------------------------------------------------------------
136: !
137: !   FormInitialGuess - Computes an initial approximation to the solution.
138: !
139: !   Input Parameters:
140: !   X    - vector
141: !
142: !   Output Parameters:
143: !   X    - vector
144: !   ierr - error code
145: !
146:       subroutine FormInitialGuess(X,ierr)
147:       use eptorsion2fmodule
148:       implicit none

150: !  Input/output variables:
151:       Vec              X
152:       PetscErrorCode   ierr

154: !  Local variables:
155:       PetscInt         i, j, k, xe, ye
156:       PetscReal      temp, val, hx, hy
157:       PetscInt         xs, ys, xm, ym
158:       PetscInt         gxm, gym, gxs, gys
159:       PetscInt         i1

161:       i1 = 1
162:       hx = 1.0/real(mx + 1)
163:       hy = 1.0/real(my + 1)

165: !  Get corner information
166:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
167:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

169: !  Compute initial guess over locally owned part of mesh
170:       xe = xs+xm
171:       ye = ys+ym
172:       do j=ys,ye-1
173:          temp = min(j+1,my-j)*hy
174:          do i=xs,xe-1
175:             k   = (j-gys)*gxm + i-gxs
176:             val = min((min(i+1,mx-i))*hx,temp)
177:             PetscCall(VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr))
178:          end do
179:       end do
180:       PetscCall(VecAssemblyBegin(X,ierr))
181:       PetscCall(VecAssemblyEnd(X,ierr))
182:       return
183:       end

185: ! ---------------------------------------------------------------------
186: !
187: !  FormFunctionGradient - Evaluates gradient G(X).
188: !
189: !  Input Parameters:
190: !  tao   - the Tao context
191: !  X     - input vector
192: !  dummy - optional user-defined context (not used here)
193: !
194: !  Output Parameters:
195: !  f     - the function value at X
196: !  G     - vector containing the newly evaluated gradient
197: !  ierr  - error code
198: !
199: !  Notes:
200: !  This routine serves as a wrapper for the lower-level routine
201: !  "ApplicationGradient", where the actual computations are
202: !  done using the standard Fortran style of treating the local
203: !  input vector data as an array over the local mesh.
204: !
205:       subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
206:       use eptorsion2fmodule
207:       implicit none

209: !  Input/output variables:
210:       Tao        tao
211:       Vec              X, G
212:       PetscReal      f
213:       PetscErrorCode   ierr
214:       PetscInt         dummy

216: !  Declarations for use with local array:

218: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
219: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
220: ! will return an array of doubles referenced by x_array offset by x_index.
221: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
222: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
223:       PetscReal      lx_v(0:1)
224:       PetscOffset      lx_i

226: !  Local variables:
227:       PetscReal      zero, p5, area, cdiv3
228:       PetscReal      val, flin, fquad,floc
229:       PetscReal      v, vb, vl, vr, vt, dvdx
230:       PetscReal      dvdy, hx, hy
231:       PetscInt         xe, ye, xsm, ysm
232:       PetscInt         xep, yep, i, j, k, ind
233:       PetscInt         xs, ys, xm, ym
234:       PetscInt         gxs, gys, gxm, gym
235:       PetscInt         i1

237:       i1 = 1
238:       0
239:       cdiv3 = param/3.0
240:       zero = 0.0
241:       p5   = 0.5
242:       hx = 1.0/real(mx + 1)
243:       hy = 1.0/real(my + 1)
244:       fquad = zero
245:       flin = zero

247: !  Initialize gradient to zero
248:       PetscCall(VecSet(G,zero,ierr))

250: !  Scatter ghost points to local vector
251:       PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
252:       PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))

254: !  Get corner information
255:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
256:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

258: !  Get pointer to vector data.
259:       PetscCall(VecGetArray(localX,lx_v,lx_i,ierr))

261: !  Set local loop dimensions
262:       xe = xs+xm
263:       ye = ys+ym
264:       if (xs .eq. 0) then
265:          xsm = xs-1
266:       else
267:          xsm = xs
268:       endif
269:       if (ys .eq. 0) then
270:          ysm = ys-1
271:       else
272:          ysm = ys
273:       endif
274:       if (xe .eq. mx) then
275:          xep = xe+1
276:       else
277:          xep = xe
278:       endif
279:       if (ye .eq. my) then
280:          yep = ye+1
281:       else
282:          yep = ye
283:       endif

285: !     Compute local gradient contributions over the lower triangular elements

287:       do j = ysm, ye-1
288:          do i = xsm, xe-1
289:             k  = (j-gys)*gxm + i-gxs
290:             v  = zero
291:             vr = zero
292:             vt = zero
293:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(lx_i+k)
294:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
295:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
296:             dvdx = (vr-v)/hx
297:             dvdy = (vt-v)/hy
298:             if (i .ne. -1 .and. j .ne. -1) then
299:                ind = k
300:                val = - dvdx/hx - dvdy/hy - cdiv3
301:                PetscCall(VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr))
302:             endif
303:             if (i .ne. mx-1 .and. j .ne. -1) then
304:                ind = k+1
305:                val =  dvdx/hx - cdiv3
306:                PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
307:             endif
308:             if (i .ne. -1 .and. j .ne. my-1) then
309:               ind = k+gxm
310:               val = dvdy/hy - cdiv3
311:               PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
312:             endif
313:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
314:             flin = flin - cdiv3 * (v+vr+vt)
315:          end do
316:       end do

318: !     Compute local gradient contributions over the upper triangular elements

320:       do j = ys, yep-1
321:          do i = xs, xep-1
322:             k  = (j-gys)*gxm + i-gxs
323:             vb = zero
324:             vl = zero
325:             v  = zero
326:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
327:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
328:             if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
329:             dvdx = (v-vl)/hx
330:             dvdy = (v-vb)/hy
331:             if (i .ne. mx .and. j .ne. 0) then
332:                ind = k-gxm
333:                val = - dvdy/hy - cdiv3
334:                PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
335:             endif
336:             if (i .ne. 0 .and. j .ne. my) then
337:                ind = k-1
338:                val =  - dvdx/hx - cdiv3
339:                PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
340:             endif
341:             if (i .ne. mx .and. j .ne. my) then
342:                ind = k
343:                val =  dvdx/hx + dvdy/hy - cdiv3
344:                PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
345:             endif
346:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
347:             flin = flin - cdiv3*(vb + vl + v)
348:          end do
349:       end do

351: !  Restore vector
352:       PetscCall(VecRestoreArray(localX,lx_v,lx_i,ierr))

354: !  Assemble gradient vector
355:       PetscCall(VecAssemblyBegin(G,ierr))
356:       PetscCall(VecAssemblyEnd(G,ierr))

358: ! Scale the gradient
359:       area = p5*hx*hy
360:       floc = area *(p5*fquad+flin)
361:       PetscCall(VecScale(G,area,ierr))

363: !  Sum function contributions from all processes
364:       PetscCallMPI(MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
365:       PetscCall(PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+16.0d0*(xep-xs)*(yep-ys),ierr))
366:       return
367:       end

369:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
370:       use eptorsion2fmodule
371:       implicit none

373:       Tao       tao
374:       Vec             X
375:       Mat             H,Hpre
376:       PetscErrorCode  ierr
377:       PetscInt        dummy

379:       PetscInt i,j,k
380:       PetscInt col(0:4),row
381:       PetscInt xs,xm,gxs,gxm
382:       PetscInt ys,ym,gys,gym
383:       PetscReal v(0:4)
384:       PetscInt i1

386:       i1 = 1

388: !     Get local grid boundaries
389:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
390:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

392:       do j=ys,ys+ym-1
393:          do i=xs,xs+xm-1
394:             row = (j-gys)*gxm + (i-gxs)

396:             k = 0
397:             if (j .gt. gys) then
398:                v(k) = -1.0
399:                col(k) = row-gxm
400:                k = k + 1
401:             endif

403:             if (i .gt. gxs) then
404:                v(k) = -1.0
405:                col(k) = row - 1
406:                k = k +1
407:             endif

409:             v(k) = 4.0
410:             col(k) = row
411:             k = k + 1

413:             if (i+1 .lt. gxs + gxm) then
414:                v(k) = -1.0
415:                col(k) = row + 1
416:                k = k + 1
417:             endif

419:             if (j+1 .lt. gys + gym) then
420:                v(k) = -1.0
421:                col(k) = row + gxm
422:                k = k + 1
423:             endif

425:             PetscCall(MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr))
426:          enddo
427:       enddo

429: !     Assemble matrix
430:       PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
431:       PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))

433: !     Tell the matrix we will never add a new nonzero location to the
434: !     matrix.  If we do it will generate an error.

436:       PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
437:       PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

439:       PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr))

441:       0
442:       return
443:       end

445:       subroutine Monitor(tao, dummy, ierr)
446:       use eptorsion2fmodule
447:       implicit none

449:       Tao tao
450:       PetscInt dummy
451:       PetscErrorCode ierr

453:       PetscInt its
454:       PetscReal f,gnorm,cnorm,xdiff
455:       TaoConvergedReason reason

457:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
458:       if (mod(its,5) .ne. 0) then
459:          PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr))
460:       endif

462:       0

464:       return
465:       end

467:       subroutine ConvergenceTest(tao, dummy, ierr)
468:       use eptorsion2fmodule
469:       implicit none

471:       Tao tao
472:       PetscInt dummy
473:       PetscErrorCode ierr

475:       PetscInt its
476:       PetscReal f,gnorm,cnorm,xdiff
477:       TaoConvergedReason reason

479:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
480:       if (its .eq. 7) then
481:        PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr))
482:       endif

484:       0

486:       return
487:       end

489: !/*TEST
490: !
491: !   build:
492: !      requires: !complex
493: !
494: !   test:
495: !      args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
496: !
497: !   test:
498: !      suffix: 2
499: !      nsize: 2
500: !      args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
501: !
502: !   test:
503: !      suffix: 3
504: !      args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
505: !TEST*/