Actual source code: eptorsion1.c
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)}
12:
13: where C is the torsion angle per unit length.
15: The multiprocessor version of this code is eptorsion2.c.
17: ---------------------------------------------------------------------- */
19: /*
20: Include "taosolver.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 "taosolver.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: TaoInitialize(); TaoFinalize();
45: Routines: TaoCreate(); TaoSetType();
46: Routines: TaoSetInitialVector();
47: Routines: TaoSetObjectiveAndGradientRoutine();
48: Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
49: Routines: TaoGetKSP(); TaoSolve();
50: Routines: TaoGetTerminationReason(); TaoDestroy();
51: Processors: 1
52: T*/
54: /*
55: User-defined application context - contains data needed by the
56: application-provided call-back routines, FormFunction(),
57: FormGradient(), and FormHessian().
58: */
60: typedef struct {
61: PetscReal param; /* nonlinearity parameter */
62: PetscInt mx, my; /* discretization in x- and y-directions */
63: PetscInt ndim; /* problem dimension */
64: Vec s, y, xvec; /* work space for computing Hessian */
65: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
66: } AppCtx;
68: /* -------- User-defined Routines --------- */
70: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
71: PetscErrorCode FormFunction(TaoSolver,Vec,PetscReal*,void*);
72: PetscErrorCode FormGradient(TaoSolver,Vec,Vec,void*);
73: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*, MatStructure *,void*);
74: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
75: PetscErrorCode HessianProduct(void*,Vec,Vec);
76: PetscErrorCode MatrixFreeHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);
77: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal *,Vec,void *);
81: PetscErrorCode main(int argc,char **argv)
82: {
84: PetscInt mx=10; /* discretization in x-direction */
85: PetscInt my=10; /* discretization in y-direction */
86: Vec x; /* solution, gradient vectors */
87: PetscBool flg; /* A return value when checking for use options */
88: TaoSolver tao; /* TaoSolver solver context */
89: Mat H; /* Hessian matrix */
90: TaoSolverTerminationReason reason;
91: KSP ksp; /* PETSc Krylov subspace solver */
92: AppCtx user; /* application context */
93: PetscMPIInt size; /* number of processes */
94: PetscReal one=1.0;
97: /* Initialize TAO,PETSc */
98: PetscInitialize(&argc,&argv,(char *)0,help);
99: TaoInitialize(&argc,&argv,(char *)0,help);
101: /* Optional: Check that only one processor is being used. */
102: MPI_Comm_size(MPI_COMM_WORLD,&size);
103: if (size >1) {
104: PetscPrintf(PETSC_COMM_SELF,"This example is intended for single processor use!\n");
105: PetscPrintf(PETSC_COMM_SELF,"Try the example eptorsion2!\n");
106: SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");
107: }
109: /* Specify default parameters for the problem, check for command-line overrides */
110: user.param = 5.0;
111: PetscOptionsGetInt(PETSC_NULL,"-my",&my,&flg);
112: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,&flg);
113: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,&flg);
116: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
117: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
118: user.ndim = mx * my; user.mx = mx; user.my = my;
120: user.hx = one/(mx+1); user.hy = one/(my+1);
123: /* Allocate vectors */
124: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
125: VecDuplicate(user.y,&user.s);
126: VecDuplicate(user.y,&x);
128: /* The TAO code begins here */
130: /* Create TAO solver and set desired solution method */
131: TaoCreate(PETSC_COMM_SELF,&tao);
132: TaoSetType(tao,"tao_lmvm");
134: /* Set solution vector with an initial guess */
135: FormInitialGuess(&user,x);
136: TaoSetInitialVector(tao,x);
138: /* Set routine for function and gradient evaluation */
139: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
141: /* From command line options, determine if using matrix-free hessian */
142: PetscOptionsHasName(PETSC_NULL,"-my_tao_mf",&flg);
143: if (flg) {
144: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,
145: user.ndim,(void*)&user,&H);
146: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat); CHKERRQ
147: (ierr);
148: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
150: TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
153: /* Set null preconditioner. Alternatively, set user-provided
154: preconditioner or explicitly form preconditioning matrix */
155: PetscOptionsSetValue("-tao_pc_type","none");
157: } else {
159: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,PETSC_NULL,&H);
160: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
162: TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
164: }
168: /* Modify the PETSc KSP structure */
169: PetscOptionsSetValue("-tao_ksp_type","cg");
171: /* Check for any TAO command line options */
172: TaoSetFromOptions(tao);
175: /* SOLVE THE APPLICATION */
176: TaoSolve(tao);
177: TaoGetKSP(tao,&ksp);
178: if (ksp) {
179: KSPView(ksp,PETSC_VIEWER_STDOUT_SELF);
180: }
182: /*
183: To View TAO solver information use
184: TaoView(tao);
185: */
187: /* Get information on termination */
188: TaoGetTerminationReason(tao,&reason);
189: if (reason <= 0){
190: PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
191: }
193: /* Free TAO data structures */
194: TaoDestroy(&tao);
196: /* Free PETSc data structures */
197: VecDestroy(&user.s);
198: VecDestroy(&user.y);
199: VecDestroy(&x);
200: MatDestroy(&H);
202: /* Finalize TAO, PETSc */
203: TaoFinalize();
204: PetscFinalize();
206: return 0;
207: }
209: /* ------------------------------------------------------------------- */
212: /*
213: FormInitialGuess - Computes an initial approximation to the solution.
215: Input Parameters:
216: . user - user-defined application context
217: . X - vector
219: Output Parameters:
220: . X - vector
221: */
222: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
223: {
224: PetscReal hx = user->hx, hy = user->hy, temp;
225: PetscReal val;
227: PetscInt i, j, k, nx = user->mx, ny = user->my;
229: /* Compute initial guess */
230: for (j=0; j<ny; j++) {
231: temp = PetscMin(j+1,ny-j)*hy;
232: for (i=0; i<nx; i++) {
233: k = nx*j + i;
234: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
235: VecSetValues(X,1,&k,&val,ADD_VALUES);
236: }
237: }
238: VecAssemblyBegin(X);
239: VecAssemblyEnd(X);
240: return 0;
241: }
242: /* ------------------------------------------------------------------- */
245: /*
246: FormFunctionGradient - Evaluates the function and corresponding gradient.
247:
248: Input Parameters:
249: tao - the TaoSolver context
250: X - the input vector
251: ptr - optional user-defined context, as set by TaoSetFunction()
253: Output Parameters:
254: f - the newly evaluated function
255: G - the newly evaluated gradient
256: */
257: PetscErrorCode FormFunctionGradient(TaoSolver tao,Vec X,PetscReal *f,Vec G,void *ptr)
258: {
260: FormFunction(tao,X,f,ptr);
261: FormGradient(tao,X,G,ptr);
262: return 0;
263: }
264: /* ------------------------------------------------------------------- */
267: /*
268: FormFunction - Evaluates the function, f(X).
270: Input Parameters:
271: . tao - the TaoSolver context
272: . X - the input vector
273: . ptr - optional user-defined context, as set by TaoSetFunction()
275: Output Parameters:
276: . f - the newly evaluated function
277: */
278: PetscErrorCode FormFunction(TaoSolver tao,Vec X,PetscReal *f,void *ptr)
279: {
280: AppCtx *user = (AppCtx *) ptr;
281: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
282: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
283: PetscReal v, cdiv3 = user->param/three;
284: PetscReal *x;
286: PetscInt nx = user->mx, ny = user->my, i, j, k;
288: /* Get pointer to vector data */
289: VecGetArray(X,&x);
291: /* Compute function contributions over the lower triangular elements */
292: for (j=-1; j<ny; j++) {
293: for (i=-1; i<nx; i++) {
294: k = nx*j + i;
295: v = zero;
296: vr = zero;
297: vt = zero;
298: if (i >= 0 && j >= 0) v = x[k];
299: if (i < nx-1 && j > -1) vr = x[k+1];
300: if (i > -1 && j < ny-1) vt = x[k+nx];
301: dvdx = (vr-v)/hx;
302: dvdy = (vt-v)/hy;
303: fquad += dvdx*dvdx + dvdy*dvdy;
304: flin -= cdiv3*(v+vr+vt);
305: }
306: }
308: /* Compute function contributions over the upper triangular elements */
309: for (j=0; j<=ny; j++) {
310: for (i=0; i<=nx; i++) {
311: k = nx*j + i;
312: vb = zero;
313: vl = zero;
314: v = zero;
315: if (i < nx && j > 0) vb = x[k-nx];
316: if (i > 0 && j < ny) vl = x[k-1];
317: if (i < nx && j < ny) v = x[k];
318: dvdx = (v-vl)/hx;
319: dvdy = (v-vb)/hy;
320: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
321: flin = flin - cdiv3*(vb+vl+v);
322: }
323: }
325: /* Restore vector */
326: VecRestoreArray(X,&x);
328: /* Scale the function */
329: area = p5*hx*hy;
330: *f = area*(p5*fquad+flin);
332: PetscLogFlops(nx*ny*24);
333: return 0;
334: }
335: /* ------------------------------------------------------------------- */
338: /*
339: FormGradient - Evaluates the gradient, G(X).
341: Input Parameters:
342: . tao - the TaoSolver context
343: . X - input vector
344: . ptr - optional user-defined context
345:
346: Output Parameters:
347: . G - vector containing the newly evaluated gradient
348: */
349: PetscErrorCode FormGradient(TaoSolver tao,Vec X,Vec G,void *ptr)
350: {
351: AppCtx *user = (AppCtx *) ptr;
352: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val, *x;
354: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
355: PetscReal hx = user->hx, hy = user->hy;
356: PetscReal vb, vl, vr, vt, dvdx, dvdy;
357: PetscReal v, cdiv3 = user->param/three;
359: /* Initialize gradient to zero */
360: VecSet(G, zero);
362: /* Get pointer to vector data */
363: VecGetArray(X,&x);
365: /* Compute gradient contributions over the lower triangular elements */
366: for (j=-1; j<ny; j++) {
367: for (i=-1; i<nx; i++) {
368: k = nx*j + i;
369: v = zero;
370: vr = zero;
371: vt = zero;
372: if (i >= 0 && j >= 0) v = x[k];
373: if (i < nx-1 && j > -1) vr = x[k+1];
374: if (i > -1 && j < ny-1) vt = x[k+nx];
375: dvdx = (vr-v)/hx;
376: dvdy = (vt-v)/hy;
377: if (i != -1 && j != -1) {
378: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
379: VecSetValues(G,1,&ind,&val,ADD_VALUES);
380: }
381: if (i != nx-1 && j != -1) {
382: ind = k+1; val = dvdx/hx - cdiv3;
383: VecSetValues(G,1,&ind,&val,ADD_VALUES);
384: }
385: if (i != -1 && j != ny-1) {
386: ind = k+nx; val = dvdy/hy - cdiv3;
387: VecSetValues(G,1,&ind,&val,ADD_VALUES);
388: }
389: }
390: }
392: /* Compute gradient contributions over the upper triangular elements */
393: for (j=0; j<=ny; j++) {
394: for (i=0; i<=nx; i++) {
395: k = nx*j + i;
396: vb = zero;
397: vl = zero;
398: v = zero;
399: if (i < nx && j > 0) vb = x[k-nx];
400: if (i > 0 && j < ny) vl = x[k-1];
401: if (i < nx && j < ny) v = x[k];
402: dvdx = (v-vl)/hx;
403: dvdy = (v-vb)/hy;
404: if (i != nx && j != 0) {
405: ind = k-nx; val = - dvdy/hy - cdiv3;
406: VecSetValues(G,1,&ind,&val,ADD_VALUES);
407: }
408: if (i != 0 && j != ny) {
409: ind = k-1; val = - dvdx/hx - cdiv3;
410: VecSetValues(G,1,&ind,&val,ADD_VALUES);
411: }
412: if (i != nx && j != ny) {
413: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
414: VecSetValues(G,1,&ind,&val,ADD_VALUES);
415: }
416: }
417: }
419: /* Restore vector */
420: VecRestoreArray(X,&x);
422: /* Assemble gradient vector */
423: VecAssemblyBegin(G);
424: VecAssemblyEnd(G);
426: /* Scale the gradient */
427: area = p5*hx*hy;
428: VecScale(G, area);
429:
430: PetscLogFlops(nx*ny*24);
431: return 0;
432: }
434: /* ------------------------------------------------------------------- */
437: /*
438: FormHessian - Forms the Hessian matrix.
440: Input Parameters:
441: . tao - the TaoSolver context
442: . X - the input vector
443: . ptr - optional user-defined context, as set by TaoSetHessian()
444:
445: Output Parameters:
446: . H - Hessian matrix
447: . PrecH - optionally different preconditioning Hessian
448: . flag - flag indicating matrix structure
450: Notes:
451: This routine is intended simply as an example of the interface
452: to a Hessian evaluation routine. Since this example compute the
453: Hessian a column at a time, it is not particularly efficient and
454: is not recommended.
455: */
456: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *HH,Mat *Hpre, MatStructure *flg, void *ptr)
457: {
458: AppCtx *user = (AppCtx *) ptr;
460: PetscInt i,j, ndim = user->ndim;
461: PetscReal *y, zero = 0.0, one = 1.0;
462: Mat H=*HH;
463: PetscBool assembled;
464: *Hpre = H;
466: /* Set location of vector */
467: user->xvec = X;
469: /* Initialize Hessian entries and work vector to zero */
470: MatAssembled(H,&assembled);
471: if (assembled){MatZeroEntries(H); }
473: VecSet(user->s, zero);
475: /* Loop over matrix columns to compute entries of the Hessian */
476: for (j=0; j<ndim; j++) {
478: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
479: VecAssemblyBegin(user->s);
480: VecAssemblyEnd(user->s);
482: HessianProduct(ptr,user->s,user->y);
484: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
485: VecAssemblyBegin(user->s);
486: VecAssemblyEnd(user->s);
488: VecGetArray(user->y,&y);
489: for (i=0; i<ndim; i++) {
490: if (y[i] != zero) {
491: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
492: }
493: }
494: VecRestoreArray(user->y,&y);
496: }
498: *flg=SAME_NONZERO_PATTERN;
500: /* Assemble matrix */
501: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
502: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
503: return 0;
504: }
507: /* ------------------------------------------------------------------- */
510: /*
511: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
512: products.
513:
514: Input Parameters:
515: . tao - the TaoSolver context
516: . X - the input vector
517: . ptr - optional user-defined context, as set by TaoSetHessian()
518:
519: Output Parameters:
520: . H - Hessian matrix
521: . PrecH - optionally different preconditioning Hessian
522: . flag - flag indicating matrix structure
523: */
524: PetscErrorCode MatrixFreeHessian(TaoSolver tao,Vec X,Mat *H,Mat *PrecH,
525: MatStructure *flag,void *ptr)
526: {
527: AppCtx *user = (AppCtx *) ptr;
529: /* Sets location of vector for use in computing matrix-vector products
530: of the form H(X)*y */
532: user->xvec = X;
533: return 0;
534: }
536: /* ------------------------------------------------------------------- */
539: /*
540: HessianProductMat - Computes the matrix-vector product
541: y = mat*svec.
543: Input Parameters:
544: . mat - input matrix
545: . svec - input vector
547: Output Parameters:
548: . y - solution vector
549: */
550: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
551: {
552: void *ptr;
554: MatShellGetContext(mat,&ptr);
555: HessianProduct(ptr,svec,y);
556:
558: return 0;
559: }
561: /* ------------------------------------------------------------------- */
564: /*
565: Hessian Product - Computes the matrix-vector product:
566: y = f''(x)*svec.
568: Input Parameters
569: . ptr - pointer to the user-defined context
570: . svec - input vector
572: Output Parameters:
573: . y - product vector
574: */
575: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
576: {
577: AppCtx *user = (AppCtx *)ptr;
578: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
579: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
581: PetscInt nx, ny, i, j, k, ind;
583: nx = user->mx;
584: ny = user->my;
585: hx = user->hx;
586: hy = user->hy;
587: hxhx = one/(hx*hx);
588: hyhy = one/(hy*hy);
590: /* Get pointers to vector data */
591: VecGetArray(user->xvec,&x);
592: VecGetArray(svec,&s);
594: /* Initialize product vector to zero */
595: VecSet(y, zero);
597: /* Compute f''(x)*s over the lower triangular elements */
598: for (j=-1; j<ny; j++) {
599: for (i=-1; i<nx; i++) {
600: k = nx*j + i;
601: v = zero;
602: vr = zero;
603: vt = zero;
604: if (i != -1 && j != -1) v = s[k];
605: if (i != nx-1 && j != -1) {
606: vr = s[k+1];
607: ind = k+1; val = hxhx*(vr-v);
608: VecSetValues(y,1,&ind,&val,ADD_VALUES);
609: }
610: if (i != -1 && j != ny-1) {
611: vt = s[k+nx];
612: ind = k+nx; val = hyhy*(vt-v);
613: VecSetValues(y,1,&ind,&val,ADD_VALUES);
614: }
615: if (i != -1 && j != -1) {
616: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
617: VecSetValues(y,1,&ind,&val,ADD_VALUES);
618: }
619: }
620: }
621:
622: /* Compute f''(x)*s over the upper triangular elements */
623: for (j=0; j<=ny; j++) {
624: for (i=0; i<=nx; i++) {
625: k = nx*j + i;
626: v = zero;
627: vl = zero;
628: vb = zero;
629: if (i != nx && j != ny) v = s[k];
630: if (i != nx && j != 0) {
631: vb = s[k-nx];
632: ind = k-nx; val = hyhy*(vb-v);
633: VecSetValues(y,1,&ind,&val,ADD_VALUES);
634: }
635: if (i != 0 && j != ny) {
636: vl = s[k-1];
637: ind = k-1; val = hxhx*(vl-v);
638: VecSetValues(y,1,&ind,&val,ADD_VALUES);
639: }
640: if (i != nx && j != ny) {
641: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
642: VecSetValues(y,1,&ind,&val,ADD_VALUES);
643: }
644: }
645: }
646: /* Restore vector data */
647: VecRestoreArray(svec,&s);
648: VecRestoreArray(user->xvec,&x);
650: /* Assemble vector */
651: VecAssemblyBegin(y);
652: VecAssemblyEnd(y);
653:
654: /* Scale resulting vector by area */
655: area = p5*hx*hy;
656: VecScale(y, area);
658: PetscLogFlops(nx*ny*18);
659:
660: return 0;
661: }