Actual source code: eptorsion1.c
petsc-3.6.1 2015-08-06
1: /* Program usage: mpirun -np 1 eptorsion1 [-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
5: Elastic-plastic torsion problem.
7: The elastic plastic torsion problem arises from the determination
8: of the stress field on an infinitely long cylindrical bar, which is
9: equivalent to the solution of the following problem:
11: min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
13: where C is the torsion angle per unit length.
15: The multiprocessor version of this code is eptorsion2.c.
17: ---------------------------------------------------------------------- */
19: /*
20: Include "petsctao.h" so that we can use TAO solvers. Note that this
21: file automatically includes files for lower-level support, such as those
22: provided by the PETSc library:
23: petsc.h - base PETSc routines petscvec.h - vectors
24: petscsys.h - sysem routines petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: */
29: #include <petsctao.h>
32: static char help[]=
33: "Demonstrates use of the TAO package to solve \n\
34: unconstrained minimization problems on a single processor. This example \n\
35: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
36: test suite.\n\
37: The command line options are:\n\
38: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
39: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
40: -par <param>, where <param> = angle of twist per unit length\n\n";
42: /*T
43: Concepts: TAO^Solving an unconstrained minimization problem
44: Routines: TaoCreate(); TaoSetType();
45: Routines: TaoSetInitialVector();
46: Routines: TaoSetObjectiveAndGradientRoutine();
47: Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
48: Routines: TaoGetKSP(); TaoSolve();
49: Routines: TaoGetConvergedReason(); TaoDestroy();
50: Processors: 1
51: T*/
53: /*
54: User-defined application context - contains data needed by the
55: application-provided call-back routines, FormFunction(),
56: FormGradient(), and FormHessian().
57: */
59: typedef struct {
60: PetscReal param; /* nonlinearity parameter */
61: PetscInt mx, my; /* discretization in x- and y-directions */
62: PetscInt ndim; /* problem dimension */
63: Vec s, y, xvec; /* work space for computing Hessian */
64: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
65: } AppCtx;
67: /* -------- User-defined Routines --------- */
69: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
70: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
71: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
72: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
73: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
74: PetscErrorCode HessianProduct(void*,Vec,Vec);
75: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
76: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
80: PetscErrorCode main(int argc,char **argv)
81: {
82: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
83: PetscInt mx=10; /* discretization in x-direction */
84: PetscInt my=10; /* discretization in y-direction */
85: Vec x; /* solution, gradient vectors */
86: PetscBool flg; /* A return value when checking for use options */
87: Tao tao; /* Tao solver context */
88: Mat H; /* Hessian matrix */
89: TaoConvergedReason reason;
90: KSP ksp; /* PETSc Krylov subspace solver */
91: AppCtx user; /* application context */
92: PetscMPIInt size; /* number of processes */
93: PetscReal one=1.0;
95: /* Initialize TAO,PETSc */
96: PetscInitialize(&argc,&argv,(char *)0,help);
98: MPI_Comm_size(MPI_COMM_WORLD,&size);
99: if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");
101: /* Specify default parameters for the problem, check for command-line overrides */
102: user.param = 5.0;
103: PetscOptionsGetInt(NULL,"-my",&my,&flg);
104: PetscOptionsGetInt(NULL,"-mx",&mx,&flg);
105: PetscOptionsGetReal(NULL,"-par",&user.param,&flg);
107: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
108: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
109: user.ndim = mx * my; user.mx = mx; user.my = my;
110: user.hx = one/(mx+1); user.hy = one/(my+1);
112: /* Allocate vectors */
113: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
114: VecDuplicate(user.y,&user.s);
115: VecDuplicate(user.y,&x);
117: /* Create TAO solver and set desired solution method */
118: TaoCreate(PETSC_COMM_SELF,&tao);
119: TaoSetType(tao,TAOLMVM);
121: /* Set solution vector with an initial guess */
122: FormInitialGuess(&user,x);
123: TaoSetInitialVector(tao,x);
125: /* Set routine for function and gradient evaluation */
126: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
128: /* From command line options, determine if using matrix-free hessian */
129: PetscOptionsHasName(NULL,"-my_tao_mf",&flg);
130: if (flg) {
131: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
132: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
133: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
135: TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
137: /* Set null preconditioner. Alternatively, set user-provided
138: preconditioner or explicitly form preconditioning matrix */
139: PetscOptionsSetValue("-pc_type","none");
140: } else {
141: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
142: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
143: TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
144: }
147: /* Check for any TAO command line options */
148: TaoSetFromOptions(tao);
150: /* SOLVE THE APPLICATION */
151: TaoSolve(tao);
152: TaoGetKSP(tao,&ksp);
153: if (ksp) {
154: KSPView(ksp,PETSC_VIEWER_STDOUT_SELF);
155: }
157: /*
158: To View TAO solver information use
159: TaoView(tao);
160: */
162: /* Get information on termination */
163: TaoGetConvergedReason(tao,&reason);
164: if (reason <= 0){
165: PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
166: }
168: TaoDestroy(&tao);
169: VecDestroy(&user.s);
170: VecDestroy(&user.y);
171: VecDestroy(&x);
172: MatDestroy(&H);
174: PetscFinalize();
175: return 0;
176: }
178: /* ------------------------------------------------------------------- */
181: /*
182: FormInitialGuess - Computes an initial approximation to the solution.
184: Input Parameters:
185: . user - user-defined application context
186: . X - vector
188: Output Parameters:
189: . X - vector
190: */
191: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
192: {
193: PetscReal hx = user->hx, hy = user->hy, temp;
194: PetscReal val;
196: PetscInt i, j, k, nx = user->mx, ny = user->my;
198: /* Compute initial guess */
199: for (j=0; j<ny; j++) {
200: temp = PetscMin(j+1,ny-j)*hy;
201: for (i=0; i<nx; i++) {
202: k = nx*j + i;
203: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
204: VecSetValues(X,1,&k,&val,ADD_VALUES);
205: }
206: }
207: VecAssemblyBegin(X);
208: VecAssemblyEnd(X);
209: return 0;
210: }
212: /* ------------------------------------------------------------------- */
215: /*
216: FormFunctionGradient - Evaluates the function and corresponding gradient.
218: Input Parameters:
219: tao - the Tao context
220: X - the input vector
221: ptr - optional user-defined context, as set by TaoSetFunction()
223: Output Parameters:
224: f - the newly evaluated function
225: G - the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
228: {
231: FormFunction(tao,X,f,ptr);
232: FormGradient(tao,X,G,ptr);
233: return 0;
234: }
236: /* ------------------------------------------------------------------- */
239: /*
240: FormFunction - Evaluates the function, f(X).
242: Input Parameters:
243: . tao - the Tao context
244: . X - the input vector
245: . ptr - optional user-defined context, as set by TaoSetFunction()
247: Output Parameters:
248: . f - the newly evaluated function
249: */
250: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
251: {
252: AppCtx *user = (AppCtx *) ptr;
253: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
254: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
255: PetscReal v, cdiv3 = user->param/three;
256: PetscReal *x;
258: PetscInt nx = user->mx, ny = user->my, i, j, k;
260: /* Get pointer to vector data */
261: VecGetArray(X,&x);
263: /* Compute function contributions over the lower triangular elements */
264: for (j=-1; j<ny; j++) {
265: for (i=-1; i<nx; i++) {
266: k = nx*j + i;
267: v = zero;
268: vr = zero;
269: vt = zero;
270: if (i >= 0 && j >= 0) v = x[k];
271: if (i < nx-1 && j > -1) vr = x[k+1];
272: if (i > -1 && j < ny-1) vt = x[k+nx];
273: dvdx = (vr-v)/hx;
274: dvdy = (vt-v)/hy;
275: fquad += dvdx*dvdx + dvdy*dvdy;
276: flin -= cdiv3*(v+vr+vt);
277: }
278: }
280: /* Compute function contributions over the upper triangular elements */
281: for (j=0; j<=ny; j++) {
282: for (i=0; i<=nx; i++) {
283: k = nx*j + i;
284: vb = zero;
285: vl = zero;
286: v = zero;
287: if (i < nx && j > 0) vb = x[k-nx];
288: if (i > 0 && j < ny) vl = x[k-1];
289: if (i < nx && j < ny) v = x[k];
290: dvdx = (v-vl)/hx;
291: dvdy = (v-vb)/hy;
292: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
293: flin = flin - cdiv3*(vb+vl+v);
294: }
295: }
297: /* Restore vector */
298: VecRestoreArray(X,&x);
300: /* Scale the function */
301: area = p5*hx*hy;
302: *f = area*(p5*fquad+flin);
304: PetscLogFlops(nx*ny*24);
305: return 0;
306: }
308: /* ------------------------------------------------------------------- */
311: /*
312: FormGradient - Evaluates the gradient, G(X).
314: Input Parameters:
315: . tao - the Tao context
316: . X - input vector
317: . ptr - optional user-defined context
319: Output Parameters:
320: . G - vector containing the newly evaluated gradient
321: */
322: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
323: {
324: AppCtx *user = (AppCtx *) ptr;
325: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val, *x;
327: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
328: PetscReal hx = user->hx, hy = user->hy;
329: PetscReal vb, vl, vr, vt, dvdx, dvdy;
330: PetscReal v, cdiv3 = user->param/three;
332: /* Initialize gradient to zero */
333: VecSet(G, zero);
335: /* Get pointer to vector data */
336: VecGetArray(X,&x);
338: /* Compute gradient contributions over the lower triangular elements */
339: for (j=-1; j<ny; j++) {
340: for (i=-1; i<nx; i++) {
341: k = nx*j + i;
342: v = zero;
343: vr = zero;
344: vt = zero;
345: if (i >= 0 && j >= 0) v = x[k];
346: if (i < nx-1 && j > -1) vr = x[k+1];
347: if (i > -1 && j < ny-1) vt = x[k+nx];
348: dvdx = (vr-v)/hx;
349: dvdy = (vt-v)/hy;
350: if (i != -1 && j != -1) {
351: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
352: VecSetValues(G,1,&ind,&val,ADD_VALUES);
353: }
354: if (i != nx-1 && j != -1) {
355: ind = k+1; val = dvdx/hx - cdiv3;
356: VecSetValues(G,1,&ind,&val,ADD_VALUES);
357: }
358: if (i != -1 && j != ny-1) {
359: ind = k+nx; val = dvdy/hy - cdiv3;
360: VecSetValues(G,1,&ind,&val,ADD_VALUES);
361: }
362: }
363: }
365: /* Compute gradient contributions over the upper triangular elements */
366: for (j=0; j<=ny; j++) {
367: for (i=0; i<=nx; i++) {
368: k = nx*j + i;
369: vb = zero;
370: vl = zero;
371: v = zero;
372: if (i < nx && j > 0) vb = x[k-nx];
373: if (i > 0 && j < ny) vl = x[k-1];
374: if (i < nx && j < ny) v = x[k];
375: dvdx = (v-vl)/hx;
376: dvdy = (v-vb)/hy;
377: if (i != nx && j != 0) {
378: ind = k-nx; val = - dvdy/hy - cdiv3;
379: VecSetValues(G,1,&ind,&val,ADD_VALUES);
380: }
381: if (i != 0 && j != ny) {
382: ind = k-1; val = - dvdx/hx - cdiv3;
383: VecSetValues(G,1,&ind,&val,ADD_VALUES);
384: }
385: if (i != nx && j != ny) {
386: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
387: VecSetValues(G,1,&ind,&val,ADD_VALUES);
388: }
389: }
390: }
391: VecRestoreArray(X,&x);
393: /* Assemble gradient vector */
394: VecAssemblyBegin(G);
395: VecAssemblyEnd(G);
397: /* Scale the gradient */
398: area = p5*hx*hy;
399: VecScale(G, area);
400: PetscLogFlops(nx*ny*24);
401: return 0;
402: }
404: /* ------------------------------------------------------------------- */
407: /*
408: FormHessian - Forms the Hessian matrix.
410: Input Parameters:
411: . tao - the Tao context
412: . X - the input vector
413: . ptr - optional user-defined context, as set by TaoSetHessian()
415: Output Parameters:
416: . H - Hessian matrix
417: . PrecH - optionally different preconditioning Hessian
418: . flag - flag indicating matrix structure
420: Notes:
421: This routine is intended simply as an example of the interface
422: to a Hessian evaluation routine. Since this example compute the
423: Hessian a column at a time, it is not particularly efficient and
424: is not recommended.
425: */
426: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
427: {
428: AppCtx *user = (AppCtx *) ptr;
430: PetscInt i,j, ndim = user->ndim;
431: PetscReal *y, zero = 0.0, one = 1.0;
432: PetscBool assembled;
434: user->xvec = X;
436: /* Initialize Hessian entries and work vector to zero */
437: MatAssembled(H,&assembled);
438: if (assembled){MatZeroEntries(H); }
440: VecSet(user->s, zero);
442: /* Loop over matrix columns to compute entries of the Hessian */
443: for (j=0; j<ndim; j++) {
444: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
445: VecAssemblyBegin(user->s);
446: VecAssemblyEnd(user->s);
448: HessianProduct(ptr,user->s,user->y);
450: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
451: VecAssemblyBegin(user->s);
452: VecAssemblyEnd(user->s);
454: VecGetArray(user->y,&y);
455: for (i=0; i<ndim; i++) {
456: if (y[i] != zero) {
457: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
458: }
459: }
460: VecRestoreArray(user->y,&y);
461: }
462: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
464: return 0;
465: }
467: /* ------------------------------------------------------------------- */
470: /*
471: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
472: products.
474: Input Parameters:
475: . tao - the Tao context
476: . X - the input vector
477: . ptr - optional user-defined context, as set by TaoSetHessian()
479: Output Parameters:
480: . H - Hessian matrix
481: . PrecH - optionally different preconditioning Hessian
482: . flag - flag indicating matrix structure
483: */
484: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
485: {
486: AppCtx *user = (AppCtx *) ptr;
488: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
489: user->xvec = X;
490: return 0;
491: }
493: /* ------------------------------------------------------------------- */
496: /*
497: HessianProductMat - Computes the matrix-vector product
498: y = mat*svec.
500: Input Parameters:
501: . mat - input matrix
502: . svec - input vector
504: Output Parameters:
505: . y - solution vector
506: */
507: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
508: {
509: void *ptr;
512: MatShellGetContext(mat,&ptr);
513: HessianProduct(ptr,svec,y);
514: return 0;
515: }
517: /* ------------------------------------------------------------------- */
520: /*
521: Hessian Product - Computes the matrix-vector product:
522: y = f''(x)*svec.
524: Input Parameters
525: . ptr - pointer to the user-defined context
526: . svec - input vector
528: Output Parameters:
529: . y - product vector
530: */
531: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
532: {
533: AppCtx *user = (AppCtx *)ptr;
534: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
535: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
537: PetscInt nx, ny, i, j, k, ind;
539: nx = user->mx;
540: ny = user->my;
541: hx = user->hx;
542: hy = user->hy;
543: hxhx = one/(hx*hx);
544: hyhy = one/(hy*hy);
546: /* Get pointers to vector data */
547: VecGetArray(user->xvec,&x);
548: VecGetArray(svec,&s);
550: /* Initialize product vector to zero */
551: VecSet(y, zero);
553: /* Compute f''(x)*s over the lower triangular elements */
554: for (j=-1; j<ny; j++) {
555: for (i=-1; i<nx; i++) {
556: k = nx*j + i;
557: v = zero;
558: vr = zero;
559: vt = zero;
560: if (i != -1 && j != -1) v = s[k];
561: if (i != nx-1 && j != -1) {
562: vr = s[k+1];
563: ind = k+1; val = hxhx*(vr-v);
564: VecSetValues(y,1,&ind,&val,ADD_VALUES);
565: }
566: if (i != -1 && j != ny-1) {
567: vt = s[k+nx];
568: ind = k+nx; val = hyhy*(vt-v);
569: VecSetValues(y,1,&ind,&val,ADD_VALUES);
570: }
571: if (i != -1 && j != -1) {
572: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
573: VecSetValues(y,1,&ind,&val,ADD_VALUES);
574: }
575: }
576: }
578: /* Compute f''(x)*s over the upper triangular elements */
579: for (j=0; j<=ny; j++) {
580: for (i=0; i<=nx; i++) {
581: k = nx*j + i;
582: v = zero;
583: vl = zero;
584: vb = zero;
585: if (i != nx && j != ny) v = s[k];
586: if (i != nx && j != 0) {
587: vb = s[k-nx];
588: ind = k-nx; val = hyhy*(vb-v);
589: VecSetValues(y,1,&ind,&val,ADD_VALUES);
590: }
591: if (i != 0 && j != ny) {
592: vl = s[k-1];
593: ind = k-1; val = hxhx*(vl-v);
594: VecSetValues(y,1,&ind,&val,ADD_VALUES);
595: }
596: if (i != nx && j != ny) {
597: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
598: VecSetValues(y,1,&ind,&val,ADD_VALUES);
599: }
600: }
601: }
602: /* Restore vector data */
603: VecRestoreArray(svec,&s);
604: VecRestoreArray(user->xvec,&x);
606: /* Assemble vector */
607: VecAssemblyBegin(y);
608: VecAssemblyEnd(y);
610: /* Scale resulting vector by area */
611: area = p5*hx*hy;
612: VecScale(y, area);
613: PetscLogFlops(nx*ny*18);
614: return 0;
615: }