Actual source code: eptorsion1.c
petsc-3.9.4 2018-09-11
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 *);
78: PetscErrorCode main(int argc,char **argv)
79: {
80: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
81: PetscInt mx=10; /* discretization in x-direction */
82: PetscInt my=10; /* discretization in y-direction */
83: Vec x; /* solution, gradient vectors */
84: PetscBool flg; /* A return value when checking for use options */
85: Tao tao; /* Tao solver context */
86: Mat H; /* Hessian matrix */
87: AppCtx user; /* application context */
88: PetscMPIInt size; /* number of processes */
89: PetscReal one=1.0;
91: /* Initialize TAO,PETSc */
92: PetscInitialize(&argc,&argv,(char *)0,help);
93: MPI_Comm_size(MPI_COMM_WORLD,&size);
94: if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");
96: /* Specify default parameters for the problem, check for command-line overrides */
97: user.param = 5.0;
98: PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
99: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
100: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
102: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
103: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
104: user.ndim = mx * my; user.mx = mx; user.my = my;
105: user.hx = one/(mx+1); user.hy = one/(my+1);
107: /* Allocate vectors */
108: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
109: VecDuplicate(user.y,&user.s);
110: VecDuplicate(user.y,&x);
112: /* Create TAO solver and set desired solution method */
113: TaoCreate(PETSC_COMM_SELF,&tao);
114: TaoSetType(tao,TAOLMVM);
116: /* Set solution vector with an initial guess */
117: FormInitialGuess(&user,x);
118: TaoSetInitialVector(tao,x);
120: /* Set routine for function and gradient evaluation */
121: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
123: /* From command line options, determine if using matrix-free hessian */
124: PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
125: if (flg) {
126: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
127: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
128: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
130: TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
132: /* Set null preconditioner. Alternatively, set user-provided
133: preconditioner or explicitly form preconditioning matrix */
134: PetscOptionsSetValue(NULL,"-pc_type","none");
135: } else {
136: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
137: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
138: TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
139: }
142: /* Check for any TAO command line options */
143: TaoSetFromOptions(tao);
145: /* SOLVE THE APPLICATION */
146: TaoSolve(tao);
148: TaoDestroy(&tao);
149: VecDestroy(&user.s);
150: VecDestroy(&user.y);
151: VecDestroy(&x);
152: MatDestroy(&H);
154: PetscFinalize();
155: return ierr;
156: }
158: /* ------------------------------------------------------------------- */
159: /*
160: FormInitialGuess - Computes an initial approximation to the solution.
162: Input Parameters:
163: . user - user-defined application context
164: . X - vector
166: Output Parameters:
167: . X - vector
168: */
169: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
170: {
171: PetscReal hx = user->hx, hy = user->hy, temp;
172: PetscReal val;
174: PetscInt i, j, k, nx = user->mx, ny = user->my;
176: /* Compute initial guess */
178: for (j=0; j<ny; j++) {
179: temp = PetscMin(j+1,ny-j)*hy;
180: for (i=0; i<nx; i++) {
181: k = nx*j + i;
182: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
183: VecSetValues(X,1,&k,&val,ADD_VALUES);
184: }
185: }
186: VecAssemblyBegin(X);
187: VecAssemblyEnd(X);
188: return(0);
189: }
191: /* ------------------------------------------------------------------- */
192: /*
193: FormFunctionGradient - Evaluates the function and corresponding gradient.
195: Input Parameters:
196: tao - the Tao context
197: X - the input vector
198: ptr - optional user-defined context, as set by TaoSetFunction()
200: Output Parameters:
201: f - the newly evaluated function
202: G - the newly evaluated gradient
203: */
204: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
205: {
209: FormFunction(tao,X,f,ptr);
210: FormGradient(tao,X,G,ptr);
211: return(0);
212: }
214: /* ------------------------------------------------------------------- */
215: /*
216: FormFunction - Evaluates the function, f(X).
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: */
226: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
227: {
228: AppCtx *user = (AppCtx *) ptr;
229: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
230: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
231: PetscReal v, cdiv3 = user->param/three;
232: const PetscScalar *x;
233: PetscErrorCode ierr;
234: PetscInt nx = user->mx, ny = user->my, i, j, k;
237: /* Get pointer to vector data */
238: VecGetArrayRead(X,&x);
240: /* Compute function contributions over the lower triangular elements */
241: for (j=-1; j<ny; j++) {
242: for (i=-1; i<nx; i++) {
243: k = nx*j + i;
244: v = zero;
245: vr = zero;
246: vt = zero;
247: if (i >= 0 && j >= 0) v = x[k];
248: if (i < nx-1 && j > -1) vr = x[k+1];
249: if (i > -1 && j < ny-1) vt = x[k+nx];
250: dvdx = (vr-v)/hx;
251: dvdy = (vt-v)/hy;
252: fquad += dvdx*dvdx + dvdy*dvdy;
253: flin -= cdiv3*(v+vr+vt);
254: }
255: }
257: /* Compute function contributions over the upper triangular elements */
258: for (j=0; j<=ny; j++) {
259: for (i=0; i<=nx; i++) {
260: k = nx*j + i;
261: vb = zero;
262: vl = zero;
263: v = zero;
264: if (i < nx && j > 0) vb = x[k-nx];
265: if (i > 0 && j < ny) vl = x[k-1];
266: if (i < nx && j < ny) v = x[k];
267: dvdx = (v-vl)/hx;
268: dvdy = (v-vb)/hy;
269: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
270: flin = flin - cdiv3*(vb+vl+v);
271: }
272: }
274: /* Restore vector */
275: VecRestoreArrayRead(X,&x);
277: /* Scale the function */
278: area = p5*hx*hy;
279: *f = area*(p5*fquad+flin);
281: PetscLogFlops(nx*ny*24);
282: return(0);
283: }
285: /* ------------------------------------------------------------------- */
286: /*
287: FormGradient - Evaluates the gradient, G(X).
289: Input Parameters:
290: . tao - the Tao context
291: . X - input vector
292: . ptr - optional user-defined context
294: Output Parameters:
295: . G - vector containing the newly evaluated gradient
296: */
297: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
298: {
299: AppCtx *user = (AppCtx *) ptr;
300: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val;
301: PetscErrorCode ierr;
302: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
303: PetscReal hx = user->hx, hy = user->hy;
304: PetscReal vb, vl, vr, vt, dvdx, dvdy;
305: PetscReal v, cdiv3 = user->param/three;
306: const PetscScalar *x;
309: /* Initialize gradient to zero */
310: VecSet(G, zero);
312: /* Get pointer to vector data */
313: VecGetArrayRead(X,&x);
315: /* Compute gradient contributions over the lower triangular elements */
316: for (j=-1; j<ny; j++) {
317: for (i=-1; i<nx; i++) {
318: k = nx*j + i;
319: v = zero;
320: vr = zero;
321: vt = zero;
322: if (i >= 0 && j >= 0) v = x[k];
323: if (i < nx-1 && j > -1) vr = x[k+1];
324: if (i > -1 && j < ny-1) vt = x[k+nx];
325: dvdx = (vr-v)/hx;
326: dvdy = (vt-v)/hy;
327: if (i != -1 && j != -1) {
328: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
329: VecSetValues(G,1,&ind,&val,ADD_VALUES);
330: }
331: if (i != nx-1 && j != -1) {
332: ind = k+1; val = dvdx/hx - cdiv3;
333: VecSetValues(G,1,&ind,&val,ADD_VALUES);
334: }
335: if (i != -1 && j != ny-1) {
336: ind = k+nx; val = dvdy/hy - cdiv3;
337: VecSetValues(G,1,&ind,&val,ADD_VALUES);
338: }
339: }
340: }
342: /* Compute gradient contributions over the upper triangular elements */
343: for (j=0; j<=ny; j++) {
344: for (i=0; i<=nx; i++) {
345: k = nx*j + i;
346: vb = zero;
347: vl = zero;
348: v = zero;
349: if (i < nx && j > 0) vb = x[k-nx];
350: if (i > 0 && j < ny) vl = x[k-1];
351: if (i < nx && j < ny) v = x[k];
352: dvdx = (v-vl)/hx;
353: dvdy = (v-vb)/hy;
354: if (i != nx && j != 0) {
355: ind = k-nx; val = - dvdy/hy - cdiv3;
356: VecSetValues(G,1,&ind,&val,ADD_VALUES);
357: }
358: if (i != 0 && j != ny) {
359: ind = k-1; val = - dvdx/hx - cdiv3;
360: VecSetValues(G,1,&ind,&val,ADD_VALUES);
361: }
362: if (i != nx && j != ny) {
363: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
364: VecSetValues(G,1,&ind,&val,ADD_VALUES);
365: }
366: }
367: }
368: VecRestoreArrayRead(X,&x);
370: /* Assemble gradient vector */
371: VecAssemblyBegin(G);
372: VecAssemblyEnd(G);
374: /* Scale the gradient */
375: area = p5*hx*hy;
376: VecScale(G, area);
377: PetscLogFlops(nx*ny*24);
378: return(0);
379: }
381: /* ------------------------------------------------------------------- */
382: /*
383: FormHessian - Forms the Hessian matrix.
385: Input Parameters:
386: . tao - the Tao context
387: . X - the input vector
388: . ptr - optional user-defined context, as set by TaoSetHessian()
390: Output Parameters:
391: . H - Hessian matrix
392: . PrecH - optionally different preconditioning Hessian
393: . flag - flag indicating matrix structure
395: Notes:
396: This routine is intended simply as an example of the interface
397: to a Hessian evaluation routine. Since this example compute the
398: Hessian a column at a time, it is not particularly efficient and
399: is not recommended.
400: */
401: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
402: {
403: AppCtx *user = (AppCtx *) ptr;
405: PetscInt i,j, ndim = user->ndim;
406: PetscReal *y, zero = 0.0, one = 1.0;
407: PetscBool assembled;
410: user->xvec = X;
412: /* Initialize Hessian entries and work vector to zero */
413: MatAssembled(H,&assembled);
414: if (assembled){MatZeroEntries(H); }
416: VecSet(user->s, zero);
418: /* Loop over matrix columns to compute entries of the Hessian */
419: for (j=0; j<ndim; j++) {
420: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
421: VecAssemblyBegin(user->s);
422: VecAssemblyEnd(user->s);
424: HessianProduct(ptr,user->s,user->y);
426: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
427: VecAssemblyBegin(user->s);
428: VecAssemblyEnd(user->s);
430: VecGetArray(user->y,&y);
431: for (i=0; i<ndim; i++) {
432: if (y[i] != zero) {
433: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
434: }
435: }
436: VecRestoreArray(user->y,&y);
437: }
438: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
439: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
440: return(0);
441: }
443: /* ------------------------------------------------------------------- */
444: /*
445: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
446: products.
448: Input Parameters:
449: . tao - the Tao context
450: . X - the input vector
451: . ptr - optional user-defined context, as set by TaoSetHessian()
453: Output Parameters:
454: . H - Hessian matrix
455: . PrecH - optionally different preconditioning Hessian
456: . flag - flag indicating matrix structure
457: */
458: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
459: {
460: AppCtx *user = (AppCtx *) ptr;
462: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
463: user->xvec = X;
464: return(0);
465: }
467: /* ------------------------------------------------------------------- */
468: /*
469: HessianProductMat - Computes the matrix-vector product
470: y = mat*svec.
472: Input Parameters:
473: . mat - input matrix
474: . svec - input vector
476: Output Parameters:
477: . y - solution vector
478: */
479: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
480: {
481: void *ptr;
485: MatShellGetContext(mat,&ptr);
486: HessianProduct(ptr,svec,y);
487: return(0);
488: }
490: /* ------------------------------------------------------------------- */
491: /*
492: Hessian Product - Computes the matrix-vector product:
493: y = f''(x)*svec.
495: Input Parameters
496: . ptr - pointer to the user-defined context
497: . svec - input vector
499: Output Parameters:
500: . y - product vector
501: */
502: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
503: {
504: AppCtx *user = (AppCtx *)ptr;
505: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
506: const PetscScalar *x, *s;
507: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
508: PetscErrorCode ierr;
509: PetscInt nx, ny, i, j, k, ind;
512: nx = user->mx;
513: ny = user->my;
514: hx = user->hx;
515: hy = user->hy;
516: hxhx = one/(hx*hx);
517: hyhy = one/(hy*hy);
519: /* Get pointers to vector data */
520: VecGetArrayRead(user->xvec,&x);
521: VecGetArrayRead(svec,&s);
523: /* Initialize product vector to zero */
524: VecSet(y, zero);
526: /* Compute f''(x)*s over the lower triangular elements */
527: for (j=-1; j<ny; j++) {
528: for (i=-1; i<nx; i++) {
529: k = nx*j + i;
530: v = zero;
531: vr = zero;
532: vt = zero;
533: if (i != -1 && j != -1) v = s[k];
534: if (i != nx-1 && j != -1) {
535: vr = s[k+1];
536: ind = k+1; val = hxhx*(vr-v);
537: VecSetValues(y,1,&ind,&val,ADD_VALUES);
538: }
539: if (i != -1 && j != ny-1) {
540: vt = s[k+nx];
541: ind = k+nx; val = hyhy*(vt-v);
542: VecSetValues(y,1,&ind,&val,ADD_VALUES);
543: }
544: if (i != -1 && j != -1) {
545: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
546: VecSetValues(y,1,&ind,&val,ADD_VALUES);
547: }
548: }
549: }
551: /* Compute f''(x)*s over the upper triangular elements */
552: for (j=0; j<=ny; j++) {
553: for (i=0; i<=nx; i++) {
554: k = nx*j + i;
555: v = zero;
556: vl = zero;
557: vb = zero;
558: if (i != nx && j != ny) v = s[k];
559: if (i != nx && j != 0) {
560: vb = s[k-nx];
561: ind = k-nx; val = hyhy*(vb-v);
562: VecSetValues(y,1,&ind,&val,ADD_VALUES);
563: }
564: if (i != 0 && j != ny) {
565: vl = s[k-1];
566: ind = k-1; val = hxhx*(vl-v);
567: VecSetValues(y,1,&ind,&val,ADD_VALUES);
568: }
569: if (i != nx && j != ny) {
570: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
571: VecSetValues(y,1,&ind,&val,ADD_VALUES);
572: }
573: }
574: }
575: /* Restore vector data */
576: VecRestoreArrayRead(svec,&s);
577: VecRestoreArrayRead(user->xvec,&x);
579: /* Assemble vector */
580: VecAssemblyBegin(y);
581: VecAssemblyEnd(y);
583: /* Scale resulting vector by area */
584: area = p5*hx*hy;
585: VecScale(y, area);
586: PetscLogFlops(nx*ny*18);
587: return(0);
588: }
591: /*TEST
593: build:
594: requires: !complex
596: test:
597: args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
599: test:
600: suffix: 2
601: args: -tao_smonitor -tao_type ntl -tao_gttol 1.e-2
603: test:
604: suffix: 3
605: args: -tao_smonitor -tao_type ntr -tao_gttol 1.e-2
607: TEST*/