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/petscdmda.h"
 28: #include "petsc/finclude/petsctao.h"
 29:   use petscdmda
 30:   use petsctao
 31:   implicit none

 33:   type(tVec) localX
 34:   type(tDM) dm
 35:   PetscReal param
 36:   PetscInt mx, my
 37: end module

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

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

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

 62: i1 = 1

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

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

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

 79: !     Set up distributed array and vectors
 80: 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))
 81: PetscCallA(DMSetFromOptions(dm, ierr))
 82: PetscCallA(DMSetUp(dm, ierr))

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

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

 92: !     The TAO code begins here

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

 98: !     Set routines for function and gradient evaluation

100: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
101: PetscCallA(TaoSetHessian(ta, H, H, ComputeHessian, 0, ierr))

103: !     Set initial guess
104: PetscCallA(FormInitialGuess(x, ierr))
105: PetscCallA(TaoSetSolution(ta, x, ierr))

107: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testmonitor', flg, ierr))
108: if (flg) then
109:   PetscCallA(TaoMonitorSet(ta, Monitor, dummy, PETSC_NULL_FUNCTION, ierr))
110: end if

112: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testconvergence', flg, ierr))
113: if (flg) then
114:   PetscCallA(TaoSetConvergenceTest(ta, ConvergenceTest, dummy, ierr))
115: end if

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

120: !     SOLVE THE APPLICATION
121: PetscCallA(TaoSolve(ta, ierr))

123: !     Free TAO data structures
124: PetscCallA(TaoDestroy(ta, ierr))

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

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

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

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

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

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

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

170: !  Compute initial guess over locally owned part of mesh
171:   xe = xs + xm
172:   ye = ys + ym
173:   do j = ys, ye - 1
174:     temp = min(j + 1, my - j)*hy
175:     do i = xs, xe - 1
176:       k = (j - gys)*gxm + i - gxs
177:       val = min((min(i + 1, mx - i))*hx, temp)
178:       PetscCall(VecSetValuesLocal(X, i1, [k], [val], ADD_VALUES, ierr))
179:     end do
180:   end do
181:   PetscCall(VecAssemblyBegin(X, ierr))
182:   PetscCall(VecAssemblyEnd(X, ierr))
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(ta, X, f, G, dummy, ierr)
206:   use eptorsion2fmodule
207:   implicit none

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

216: !  Declarations for use with local array:

218:   PetscReal, pointer :: lx_v(:)

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

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

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

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

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

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

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

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

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

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

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

345: !  Restore vector
346:   PetscCall(VecRestoreArrayRead(localX, lx_v, ierr))

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

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

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

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

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

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

379:   i1 = 1

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

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

389:       k = 0
390:       if (j > gys) then
391:         v(k) = -1.0
392:         col(k) = row - gxm
393:         k = k + 1
394:       end if

396:       if (i > gxs) then
397:         v(k) = -1.0
398:         col(k) = row - 1
399:         k = k + 1
400:       end if

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

406:       if (i + 1 < gxs + gxm) then
407:         v(k) = -1.0
408:         col(k) = row + 1
409:         k = k + 1
410:       end if

412:       if (j + 1 < gys + gym) then
413:         v(k) = -1.0
414:         col(k) = row + gxm
415:         k = k + 1
416:       end if

418:       PetscCall(MatSetValuesLocal(H, i1, [row], k, col, v, INSERT_VALUES, ierr))
419:     end do
420:   end do

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

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

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

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

434:   ierr = 0
435: end

437: subroutine Monitor(ta, dummy, ierr)
438:   use eptorsion2fmodule
439:   implicit none

441:   type(tTao) ta
442:   PetscInt dummy
443:   PetscErrorCode ierr

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

449:   PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
450:   if (mod(its, 5) /= 0) then
451:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, 'iteration multiple of 5\n', ierr))
452:   end if

454:   ierr = 0

456: end

458: subroutine ConvergenceTest(ta, dummy, ierr)
459:   use eptorsion2fmodule
460:   implicit none

462:   type(tTao) ta
463:   PetscInt dummy
464:   PetscErrorCode ierr

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

470:   PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
471:   if (its == 7) then
472:     PetscCall(TaoSetConvergedReason(ta, TAO_DIVERGED_MAXITS, ierr))
473:   end if

475:   ierr = 0

477: end

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