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:       type(tVec)       localX
 33:       type(tDM)        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:       type(tVec)       x              ! solution vector
 48:       type(tMat)       H              ! hessian matrix
 49:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 50:       type(tTao)       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_ARRAY,PETSC_NULL_INTEGER_ARRAY,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(TaoMonitorSet(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:       end

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

208: !  Input/output variables:
209:       type(tTao)       tao
210:       type(tVec)       X, G
211:       PetscReal        f
212:       PetscErrorCode   ierr
213:       PetscInt         dummy

215: !  Declarations for use with local array:

217:       PetscReal, pointer :: lx_v(:)

219: !  Local variables:
220:       PetscReal      zero, p5, area, cdiv3
221:       PetscReal      val, flin, fquad,floc
222:       PetscReal      v, vb, vl, vr, vt, dvdx
223:       PetscReal      dvdy, hx, hy
224:       PetscInt       xe, ye, xsm, ysm
225:       PetscInt       xep, yep, i, j, k, ind
226:       PetscInt       xs, ys, xm, ym
227:       PetscInt       gxs, gys, gxm, gym
228:       PetscInt       i1

230:       i1 = 1
231:       ierr  = 0
232:       cdiv3 = param/3.0
233:       zero = 0.0
234:       p5   = 0.5
235:       hx = 1.0/real(mx + 1)
236:       hy = 1.0/real(my + 1)
237:       fquad = zero
238:       flin = zero

240: !  Initialize gradient to zero
241:       PetscCall(VecSet(G,zero,ierr))

243: !  Scatter ghost points to local vector
244:       PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
245:       PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))

247: !  Get corner information
248:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
249:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

251: !  Get pointer to vector data.
252:       PetscCall(VecGetArrayReadF90(localX,lx_v,ierr))

254: !  Set local loop dimensions
255:       xe = xs+xm
256:       ye = ys+ym
257:       if (xs .eq. 0) then
258:          xsm = xs-1
259:       else
260:          xsm = xs
261:       endif
262:       if (ys .eq. 0) then
263:          ysm = ys-1
264:       else
265:          ysm = ys
266:       endif
267:       if (xe .eq. mx) then
268:          xep = xe+1
269:       else
270:          xep = xe
271:       endif
272:       if (ye .eq. my) then
273:          yep = ye+1
274:       else
275:          yep = ye
276:       endif

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

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

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

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

344: !  Restore vector
345:       PetscCall(VecRestoreArrayReadF90(localX,lx_v,ierr))

347: !  Assemble gradient vector
348:       PetscCall(VecAssemblyBegin(G,ierr))
349:       PetscCall(VecAssemblyEnd(G,ierr))

351: ! Scale the gradient
352:       area = p5*hx*hy
353:       floc = area *(p5*fquad+flin)
354:       PetscCall(VecScale(G,area,ierr))

356: !  Sum function contributions from all processes
357:       PetscCallMPI(MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
358:       PetscCall(PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+16.0d0*(xep-xs)*(yep-ys),ierr))
359:       end

361:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
362:       use eptorsion2fmodule
363:       implicit none

365:       type(tTao)      tao
366:       type(tVec)      X
367:       type(tMat)      H,Hpre
368:       PetscErrorCode  ierr
369:       PetscInt        dummy

371:       PetscInt i,j,k
372:       PetscInt col(0:4),row
373:       PetscInt xs,xm,gxs,gxm
374:       PetscInt ys,ym,gys,gym
375:       PetscReal v(0:4)
376:       PetscInt i1

378:       i1 = 1

380: !     Get local grid boundaries
381:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
382:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

384:       do j=ys,ys+ym-1
385:          do i=xs,xs+xm-1
386:             row = (j-gys)*gxm + (i-gxs)

388:             k = 0
389:             if (j .gt. gys) then
390:                v(k) = -1.0
391:                col(k) = row-gxm
392:                k = k + 1
393:             endif

395:             if (i .gt. gxs) then
396:                v(k) = -1.0
397:                col(k) = row - 1
398:                k = k +1
399:             endif

401:             v(k) = 4.0
402:             col(k) = row
403:             k = k + 1

405:             if (i+1 .lt. gxs + gxm) then
406:                v(k) = -1.0
407:                col(k) = row + 1
408:                k = k + 1
409:             endif

411:             if (j+1 .lt. gys + gym) then
412:                v(k) = -1.0
413:                col(k) = row + gxm
414:                k = k + 1
415:             endif

417:             PetscCall(MatSetValuesLocal(H,i1,[row],k,col,v,INSERT_VALUES,ierr))
418:          enddo
419:       enddo

421: !     Assemble matrix
422:       PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
423:       PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))

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

428:       PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
429:       PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

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

433:       ierr = 0
434:       end

436:       subroutine Monitor(tao, dummy, ierr)
437:       use eptorsion2fmodule
438:       implicit none

440:       type(tTao)        tao
441:       PetscInt          dummy
442:       PetscErrorCode    ierr

444:       PetscInt           its
445:       PetscReal          f,gnorm,cnorm,xdiff
446:       TaoConvergedReason reason

448:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
449:       if (mod(its,5) .ne. 0) then
450:          PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr))
451:       endif

453:       ierr = 0

455:       end

457:       subroutine ConvergenceTest(tao, dummy, ierr)
458:       use eptorsion2fmodule
459:       implicit none

461:       type(tTao)          tao
462:       PetscInt           dummy
463:       PetscErrorCode     ierr

465:       PetscInt           its
466:       PetscReal          f,gnorm,cnorm,xdiff
467:       TaoConvergedReason reason

469:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
470:       if (its .eq. 7) then
471:        PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr))
472:       endif

474:       ierr = 0

476:       end

478: !/*TEST
479: !
480: !   build:
481: !      requires: !complex
482: !
483: !   test:
484: !      args: -tao_monitor_short -tao_type nls -tao_gttol 1.e-2
485: !
486: !   test:
487: !      suffix: 2
488: !      nsize: 2
489: !      args: -tao_monitor_short -tao_type lmvm -tao_gttol 1.e-2
490: !
491: !   test:
492: !      suffix: 3
493: !      args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
494: !TEST*/