Actual source code: eptorsion1.c
petsc-3.7.7 2017-09-25
1: /* Program usage: mpiexec -n 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: 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: AppCtx user; /* application context */
90: PetscMPIInt size; /* number of processes */
91: PetscReal one=1.0;
93: /* Initialize TAO,PETSc */
94: PetscInitialize(&argc,&argv,(char *)0,help);
96: MPI_Comm_size(MPI_COMM_WORLD,&size);
97: if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");
99: /* Specify default parameters for the problem, check for command-line overrides */
100: user.param = 5.0;
101: PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
102: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
103: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
105: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
106: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
107: user.ndim = mx * my; user.mx = mx; user.my = my;
108: user.hx = one/(mx+1); user.hy = one/(my+1);
110: /* Allocate vectors */
111: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
112: VecDuplicate(user.y,&user.s);
113: VecDuplicate(user.y,&x);
115: /* Create TAO solver and set desired solution method */
116: TaoCreate(PETSC_COMM_SELF,&tao);
117: TaoSetType(tao,TAOLMVM);
119: /* Set solution vector with an initial guess */
120: FormInitialGuess(&user,x);
121: TaoSetInitialVector(tao,x);
123: /* Set routine for function and gradient evaluation */
124: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
126: /* From command line options, determine if using matrix-free hessian */
127: PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
128: if (flg) {
129: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
130: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
131: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
133: TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
135: /* Set null preconditioner. Alternatively, set user-provided
136: preconditioner or explicitly form preconditioning matrix */
137: PetscOptionsSetValue(NULL,"-pc_type","none");
138: } else {
139: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
140: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
141: TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
142: }
145: /* Check for any TAO command line options */
146: TaoSetFromOptions(tao);
148: /* SOLVE THE APPLICATION */
149: TaoSolve(tao);
151: TaoDestroy(&tao);
152: VecDestroy(&user.s);
153: VecDestroy(&user.y);
154: VecDestroy(&x);
155: MatDestroy(&H);
157: PetscFinalize();
158: return 0;
159: }
161: /* ------------------------------------------------------------------- */
164: /*
165: FormInitialGuess - Computes an initial approximation to the solution.
167: Input Parameters:
168: . user - user-defined application context
169: . X - vector
171: Output Parameters:
172: . X - vector
173: */
174: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
175: {
176: PetscReal hx = user->hx, hy = user->hy, temp;
177: PetscReal val;
179: PetscInt i, j, k, nx = user->mx, ny = user->my;
181: /* Compute initial guess */
182: for (j=0; j<ny; j++) {
183: temp = PetscMin(j+1,ny-j)*hy;
184: for (i=0; i<nx; i++) {
185: k = nx*j + i;
186: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
187: VecSetValues(X,1,&k,&val,ADD_VALUES);
188: }
189: }
190: VecAssemblyBegin(X);
191: VecAssemblyEnd(X);
192: return 0;
193: }
195: /* ------------------------------------------------------------------- */
198: /*
199: FormFunctionGradient - Evaluates the function and corresponding gradient.
201: Input Parameters:
202: tao - the Tao context
203: X - the input vector
204: ptr - optional user-defined context, as set by TaoSetFunction()
206: Output Parameters:
207: f - the newly evaluated function
208: G - the newly evaluated gradient
209: */
210: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
211: {
214: FormFunction(tao,X,f,ptr);
215: FormGradient(tao,X,G,ptr);
216: return 0;
217: }
219: /* ------------------------------------------------------------------- */
222: /*
223: FormFunction - Evaluates the function, f(X).
225: Input Parameters:
226: . tao - the Tao context
227: . X - the input vector
228: . ptr - optional user-defined context, as set by TaoSetFunction()
230: Output Parameters:
231: . f - the newly evaluated function
232: */
233: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
234: {
235: AppCtx *user = (AppCtx *) ptr;
236: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
237: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
238: PetscReal v, cdiv3 = user->param/three;
239: PetscReal *x;
241: PetscInt nx = user->mx, ny = user->my, i, j, k;
243: /* Get pointer to vector data */
244: VecGetArray(X,&x);
246: /* Compute function contributions over the lower triangular elements */
247: for (j=-1; j<ny; j++) {
248: for (i=-1; i<nx; i++) {
249: k = nx*j + i;
250: v = zero;
251: vr = zero;
252: vt = zero;
253: if (i >= 0 && j >= 0) v = x[k];
254: if (i < nx-1 && j > -1) vr = x[k+1];
255: if (i > -1 && j < ny-1) vt = x[k+nx];
256: dvdx = (vr-v)/hx;
257: dvdy = (vt-v)/hy;
258: fquad += dvdx*dvdx + dvdy*dvdy;
259: flin -= cdiv3*(v+vr+vt);
260: }
261: }
263: /* Compute function contributions over the upper triangular elements */
264: for (j=0; j<=ny; j++) {
265: for (i=0; i<=nx; i++) {
266: k = nx*j + i;
267: vb = zero;
268: vl = zero;
269: v = zero;
270: if (i < nx && j > 0) vb = x[k-nx];
271: if (i > 0 && j < ny) vl = x[k-1];
272: if (i < nx && j < ny) v = x[k];
273: dvdx = (v-vl)/hx;
274: dvdy = (v-vb)/hy;
275: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
276: flin = flin - cdiv3*(vb+vl+v);
277: }
278: }
280: /* Restore vector */
281: VecRestoreArray(X,&x);
283: /* Scale the function */
284: area = p5*hx*hy;
285: *f = area*(p5*fquad+flin);
287: PetscLogFlops(nx*ny*24);
288: return 0;
289: }
291: /* ------------------------------------------------------------------- */
294: /*
295: FormGradient - Evaluates the gradient, G(X).
297: Input Parameters:
298: . tao - the Tao context
299: . X - input vector
300: . ptr - optional user-defined context
302: Output Parameters:
303: . G - vector containing the newly evaluated gradient
304: */
305: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
306: {
307: AppCtx *user = (AppCtx *) ptr;
308: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val, *x;
310: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
311: PetscReal hx = user->hx, hy = user->hy;
312: PetscReal vb, vl, vr, vt, dvdx, dvdy;
313: PetscReal v, cdiv3 = user->param/three;
315: /* Initialize gradient to zero */
316: VecSet(G, zero);
318: /* Get pointer to vector data */
319: VecGetArray(X,&x);
321: /* Compute gradient contributions over the lower triangular elements */
322: for (j=-1; j<ny; j++) {
323: for (i=-1; i<nx; i++) {
324: k = nx*j + i;
325: v = zero;
326: vr = zero;
327: vt = zero;
328: if (i >= 0 && j >= 0) v = x[k];
329: if (i < nx-1 && j > -1) vr = x[k+1];
330: if (i > -1 && j < ny-1) vt = x[k+nx];
331: dvdx = (vr-v)/hx;
332: dvdy = (vt-v)/hy;
333: if (i != -1 && j != -1) {
334: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
335: VecSetValues(G,1,&ind,&val,ADD_VALUES);
336: }
337: if (i != nx-1 && j != -1) {
338: ind = k+1; val = dvdx/hx - cdiv3;
339: VecSetValues(G,1,&ind,&val,ADD_VALUES);
340: }
341: if (i != -1 && j != ny-1) {
342: ind = k+nx; val = dvdy/hy - cdiv3;
343: VecSetValues(G,1,&ind,&val,ADD_VALUES);
344: }
345: }
346: }
348: /* Compute gradient contributions over the upper triangular elements */
349: for (j=0; j<=ny; j++) {
350: for (i=0; i<=nx; i++) {
351: k = nx*j + i;
352: vb = zero;
353: vl = zero;
354: v = zero;
355: if (i < nx && j > 0) vb = x[k-nx];
356: if (i > 0 && j < ny) vl = x[k-1];
357: if (i < nx && j < ny) v = x[k];
358: dvdx = (v-vl)/hx;
359: dvdy = (v-vb)/hy;
360: if (i != nx && j != 0) {
361: ind = k-nx; val = - dvdy/hy - cdiv3;
362: VecSetValues(G,1,&ind,&val,ADD_VALUES);
363: }
364: if (i != 0 && j != ny) {
365: ind = k-1; val = - dvdx/hx - cdiv3;
366: VecSetValues(G,1,&ind,&val,ADD_VALUES);
367: }
368: if (i != nx && j != ny) {
369: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
370: VecSetValues(G,1,&ind,&val,ADD_VALUES);
371: }
372: }
373: }
374: VecRestoreArray(X,&x);
376: /* Assemble gradient vector */
377: VecAssemblyBegin(G);
378: VecAssemblyEnd(G);
380: /* Scale the gradient */
381: area = p5*hx*hy;
382: VecScale(G, area);
383: PetscLogFlops(nx*ny*24);
384: return 0;
385: }
387: /* ------------------------------------------------------------------- */
390: /*
391: FormHessian - Forms the Hessian matrix.
393: Input Parameters:
394: . tao - the Tao context
395: . X - the input vector
396: . ptr - optional user-defined context, as set by TaoSetHessian()
398: Output Parameters:
399: . H - Hessian matrix
400: . PrecH - optionally different preconditioning Hessian
401: . flag - flag indicating matrix structure
403: Notes:
404: This routine is intended simply as an example of the interface
405: to a Hessian evaluation routine. Since this example compute the
406: Hessian a column at a time, it is not particularly efficient and
407: is not recommended.
408: */
409: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
410: {
411: AppCtx *user = (AppCtx *) ptr;
413: PetscInt i,j, ndim = user->ndim;
414: PetscReal *y, zero = 0.0, one = 1.0;
415: PetscBool assembled;
417: user->xvec = X;
419: /* Initialize Hessian entries and work vector to zero */
420: MatAssembled(H,&assembled);
421: if (assembled){MatZeroEntries(H); }
423: VecSet(user->s, zero);
425: /* Loop over matrix columns to compute entries of the Hessian */
426: for (j=0; j<ndim; j++) {
427: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
428: VecAssemblyBegin(user->s);
429: VecAssemblyEnd(user->s);
431: HessianProduct(ptr,user->s,user->y);
433: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
434: VecAssemblyBegin(user->s);
435: VecAssemblyEnd(user->s);
437: VecGetArray(user->y,&y);
438: for (i=0; i<ndim; i++) {
439: if (y[i] != zero) {
440: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
441: }
442: }
443: VecRestoreArray(user->y,&y);
444: }
445: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
446: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
447: return 0;
448: }
450: /* ------------------------------------------------------------------- */
453: /*
454: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
455: products.
457: Input Parameters:
458: . tao - the Tao context
459: . X - the input vector
460: . ptr - optional user-defined context, as set by TaoSetHessian()
462: Output Parameters:
463: . H - Hessian matrix
464: . PrecH - optionally different preconditioning Hessian
465: . flag - flag indicating matrix structure
466: */
467: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
468: {
469: AppCtx *user = (AppCtx *) ptr;
471: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
472: user->xvec = X;
473: return 0;
474: }
476: /* ------------------------------------------------------------------- */
479: /*
480: HessianProductMat - Computes the matrix-vector product
481: y = mat*svec.
483: Input Parameters:
484: . mat - input matrix
485: . svec - input vector
487: Output Parameters:
488: . y - solution vector
489: */
490: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
491: {
492: void *ptr;
495: MatShellGetContext(mat,&ptr);
496: HessianProduct(ptr,svec,y);
497: return 0;
498: }
500: /* ------------------------------------------------------------------- */
503: /*
504: Hessian Product - Computes the matrix-vector product:
505: y = f''(x)*svec.
507: Input Parameters
508: . ptr - pointer to the user-defined context
509: . svec - input vector
511: Output Parameters:
512: . y - product vector
513: */
514: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
515: {
516: AppCtx *user = (AppCtx *)ptr;
517: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
518: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
520: PetscInt nx, ny, i, j, k, ind;
522: nx = user->mx;
523: ny = user->my;
524: hx = user->hx;
525: hy = user->hy;
526: hxhx = one/(hx*hx);
527: hyhy = one/(hy*hy);
529: /* Get pointers to vector data */
530: VecGetArray(user->xvec,&x);
531: VecGetArray(svec,&s);
533: /* Initialize product vector to zero */
534: VecSet(y, zero);
536: /* Compute f''(x)*s over the lower triangular elements */
537: for (j=-1; j<ny; j++) {
538: for (i=-1; i<nx; i++) {
539: k = nx*j + i;
540: v = zero;
541: vr = zero;
542: vt = zero;
543: if (i != -1 && j != -1) v = s[k];
544: if (i != nx-1 && j != -1) {
545: vr = s[k+1];
546: ind = k+1; val = hxhx*(vr-v);
547: VecSetValues(y,1,&ind,&val,ADD_VALUES);
548: }
549: if (i != -1 && j != ny-1) {
550: vt = s[k+nx];
551: ind = k+nx; val = hyhy*(vt-v);
552: VecSetValues(y,1,&ind,&val,ADD_VALUES);
553: }
554: if (i != -1 && j != -1) {
555: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
556: VecSetValues(y,1,&ind,&val,ADD_VALUES);
557: }
558: }
559: }
561: /* Compute f''(x)*s over the upper triangular elements */
562: for (j=0; j<=ny; j++) {
563: for (i=0; i<=nx; i++) {
564: k = nx*j + i;
565: v = zero;
566: vl = zero;
567: vb = zero;
568: if (i != nx && j != ny) v = s[k];
569: if (i != nx && j != 0) {
570: vb = s[k-nx];
571: ind = k-nx; val = hyhy*(vb-v);
572: VecSetValues(y,1,&ind,&val,ADD_VALUES);
573: }
574: if (i != 0 && j != ny) {
575: vl = s[k-1];
576: ind = k-1; val = hxhx*(vl-v);
577: VecSetValues(y,1,&ind,&val,ADD_VALUES);
578: }
579: if (i != nx && j != ny) {
580: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
581: VecSetValues(y,1,&ind,&val,ADD_VALUES);
582: }
583: }
584: }
585: /* Restore vector data */
586: VecRestoreArray(svec,&s);
587: VecRestoreArray(user->xvec,&x);
589: /* Assemble vector */
590: VecAssemblyBegin(y);
591: VecAssemblyEnd(y);
593: /* Scale resulting vector by area */
594: area = p5*hx*hy;
595: VecScale(y, area);
596: PetscLogFlops(nx*ny*18);
597: return 0;
598: }