Actual source code: eptorsion1.c
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 - system 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>
31: static char help[]=
32: "Demonstrates use of the TAO package to solve \n\
33: unconstrained minimization problems on a single processor. This example \n\
34: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
35: test suite.\n\
36: The command line options are:\n\
37: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
38: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
39: -par <param>, where <param> = angle of twist per unit length\n\n";
41: /*T
42: Concepts: TAO^Solving an unconstrained minimization problem
43: Routines: TaoCreate(); TaoSetType();
44: Routines: TaoSetSolution();
45: Routines: TaoSetObjectiveAndGradient();
46: Routines: TaoSetHessian(); TaoSetFromOptions();
47: Routines: TaoGetKSP(); TaoSolve();
48: Routines: TaoDestroy();
49: Processors: 1
50: T*/
52: /*
53: User-defined application context - contains data needed by the
54: application-provided call-back routines, FormFunction(),
55: FormGradient(), and FormHessian().
56: */
58: typedef struct {
59: PetscReal param; /* nonlinearity parameter */
60: PetscInt mx, my; /* discretization in x- and y-directions */
61: PetscInt ndim; /* problem dimension */
62: Vec s, y, xvec; /* work space for computing Hessian */
63: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
64: } AppCtx;
66: /* -------- User-defined Routines --------- */
68: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
69: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
70: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
71: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
72: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
73: PetscErrorCode HessianProduct(void*,Vec,Vec);
74: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
75: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
77: PetscErrorCode main(int argc,char **argv)
78: {
79: PetscInt mx=10; /* discretization in x-direction */
80: PetscInt my=10; /* discretization in y-direction */
81: Vec x; /* solution, gradient vectors */
82: PetscBool flg; /* A return value when checking for use options */
83: Tao tao; /* Tao solver context */
84: Mat H; /* Hessian matrix */
85: AppCtx user; /* application context */
86: PetscMPIInt size; /* number of processes */
87: PetscReal one=1.0;
89: PetscBool test_lmvm = PETSC_FALSE;
90: KSP ksp;
91: PC pc;
92: Mat M;
93: Vec in, out, out2;
94: PetscReal mult_solve_dist;
96: /* Initialize TAO,PETSc */
97: PetscInitialize(&argc,&argv,(char *)0,help);
98: MPI_Comm_size(MPI_COMM_WORLD,&size);
101: /* Specify default parameters for the problem, check for command-line overrides */
102: user.param = 5.0;
103: PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
104: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
105: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
106: PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);
108: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
109: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
110: user.ndim = mx * my; user.mx = mx; user.my = my;
111: user.hx = one/(mx+1); user.hy = one/(my+1);
113: /* Allocate vectors */
114: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
115: VecDuplicate(user.y,&user.s);
116: VecDuplicate(user.y,&x);
118: /* Create TAO solver and set desired solution method */
119: TaoCreate(PETSC_COMM_SELF,&tao);
120: TaoSetType(tao,TAOLMVM);
122: /* Set solution vector with an initial guess */
123: FormInitialGuess(&user,x);
124: TaoSetSolution(tao,x);
126: /* Set routine for function and gradient evaluation */
127: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
129: /* From command line options, determine if using matrix-free hessian */
130: PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
131: if (flg) {
132: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
133: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
134: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
136: TaoSetHessian(tao,H,H,MatrixFreeHessian,(void *)&user);
137: } else {
138: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
139: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
140: TaoSetHessian(tao,H,H,FormHessian,(void *)&user);
141: }
143: /* Test the LMVM matrix */
144: if (test_lmvm) {
145: PetscOptionsSetValue(NULL, "-tao_type", "bntr");
146: PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
147: }
149: /* Check for any TAO command line options */
150: TaoSetFromOptions(tao);
152: /* SOLVE THE APPLICATION */
153: TaoSolve(tao);
155: /* Test the LMVM matrix */
156: if (test_lmvm) {
157: TaoGetKSP(tao, &ksp);
158: KSPGetPC(ksp, &pc);
159: PCLMVMGetMatLMVM(pc, &M);
160: VecDuplicate(x, &in);
161: VecDuplicate(x, &out);
162: VecDuplicate(x, &out2);
163: VecSet(in, 5.0);
164: MatMult(M, in, out);
165: MatSolve(M, out, out2);
166: VecAXPY(out2, -1.0, in);
167: VecNorm(out2, NORM_2, &mult_solve_dist);
168: PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist);
169: VecDestroy(&in);
170: VecDestroy(&out);
171: VecDestroy(&out2);
172: }
174: TaoDestroy(&tao);
175: VecDestroy(&user.s);
176: VecDestroy(&user.y);
177: VecDestroy(&x);
178: MatDestroy(&H);
180: PetscFinalize();
181: return 0;
182: }
184: /* ------------------------------------------------------------------- */
185: /*
186: FormInitialGuess - Computes an initial approximation to the solution.
188: Input Parameters:
189: . user - user-defined application context
190: . X - vector
192: Output Parameters:
193: . X - vector
194: */
195: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
196: {
197: PetscReal hx = user->hx, hy = user->hy, temp;
198: PetscReal val;
199: PetscInt i, j, k, nx = user->mx, ny = user->my;
201: /* Compute initial guess */
203: for (j=0; j<ny; j++) {
204: temp = PetscMin(j+1,ny-j)*hy;
205: for (i=0; i<nx; i++) {
206: k = nx*j + i;
207: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
208: VecSetValues(X,1,&k,&val,ADD_VALUES);
209: }
210: }
211: VecAssemblyBegin(X);
212: VecAssemblyEnd(X);
213: return 0;
214: }
216: /* ------------------------------------------------------------------- */
217: /*
218: FormFunctionGradient - Evaluates the function and corresponding gradient.
220: Input Parameters:
221: tao - the Tao context
222: X - the input vector
223: ptr - optional user-defined context, as set by TaoSetFunction()
225: Output Parameters:
226: f - the newly evaluated function
227: G - the newly evaluated gradient
228: */
229: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
230: {
232: FormFunction(tao,X,f,ptr);
233: FormGradient(tao,X,G,ptr);
234: return 0;
235: }
237: /* ------------------------------------------------------------------- */
238: /*
239: FormFunction - Evaluates the function, f(X).
241: Input Parameters:
242: . tao - the Tao context
243: . X - the input vector
244: . ptr - optional user-defined context, as set by TaoSetFunction()
246: Output Parameters:
247: . f - the newly evaluated function
248: */
249: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
250: {
251: AppCtx *user = (AppCtx *) ptr;
252: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
253: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
254: PetscReal v, cdiv3 = user->param/three;
255: const PetscScalar *x;
256: PetscInt nx = user->mx, ny = user->my, i, j, k;
259: /* Get pointer to vector data */
260: VecGetArrayRead(X,&x);
262: /* Compute function contributions over the lower triangular elements */
263: for (j=-1; j<ny; j++) {
264: for (i=-1; i<nx; i++) {
265: k = nx*j + i;
266: v = zero;
267: vr = zero;
268: vt = zero;
269: if (i >= 0 && j >= 0) v = x[k];
270: if (i < nx-1 && j > -1) vr = x[k+1];
271: if (i > -1 && j < ny-1) vt = x[k+nx];
272: dvdx = (vr-v)/hx;
273: dvdy = (vt-v)/hy;
274: fquad += dvdx*dvdx + dvdy*dvdy;
275: flin -= cdiv3*(v+vr+vt);
276: }
277: }
279: /* Compute function contributions over the upper triangular elements */
280: for (j=0; j<=ny; j++) {
281: for (i=0; i<=nx; i++) {
282: k = nx*j + i;
283: vb = zero;
284: vl = zero;
285: v = zero;
286: if (i < nx && j > 0) vb = x[k-nx];
287: if (i > 0 && j < ny) vl = x[k-1];
288: if (i < nx && j < ny) v = x[k];
289: dvdx = (v-vl)/hx;
290: dvdy = (v-vb)/hy;
291: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
292: flin = flin - cdiv3*(vb+vl+v);
293: }
294: }
296: /* Restore vector */
297: VecRestoreArrayRead(X,&x);
299: /* Scale the function */
300: area = p5*hx*hy;
301: *f = area*(p5*fquad+flin);
303: PetscLogFlops(24.0*nx*ny);
304: return 0;
305: }
307: /* ------------------------------------------------------------------- */
308: /*
309: FormGradient - Evaluates the gradient, G(X).
311: Input Parameters:
312: . tao - the Tao context
313: . X - input vector
314: . ptr - optional user-defined context
316: Output Parameters:
317: . G - vector containing the newly evaluated gradient
318: */
319: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
320: {
321: AppCtx *user = (AppCtx *) ptr;
322: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val;
323: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
324: PetscReal hx = user->hx, hy = user->hy;
325: PetscReal vb, vl, vr, vt, dvdx, dvdy;
326: PetscReal v, cdiv3 = user->param/three;
327: const PetscScalar *x;
330: /* Initialize gradient to zero */
331: VecSet(G, zero);
333: /* Get pointer to vector data */
334: VecGetArrayRead(X,&x);
336: /* Compute gradient contributions over the lower triangular elements */
337: for (j=-1; j<ny; j++) {
338: for (i=-1; i<nx; i++) {
339: k = nx*j + i;
340: v = zero;
341: vr = zero;
342: vt = zero;
343: if (i >= 0 && j >= 0) v = x[k];
344: if (i < nx-1 && j > -1) vr = x[k+1];
345: if (i > -1 && j < ny-1) vt = x[k+nx];
346: dvdx = (vr-v)/hx;
347: dvdy = (vt-v)/hy;
348: if (i != -1 && j != -1) {
349: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
350: VecSetValues(G,1,&ind,&val,ADD_VALUES);
351: }
352: if (i != nx-1 && j != -1) {
353: ind = k+1; val = dvdx/hx - cdiv3;
354: VecSetValues(G,1,&ind,&val,ADD_VALUES);
355: }
356: if (i != -1 && j != ny-1) {
357: ind = k+nx; val = dvdy/hy - cdiv3;
358: VecSetValues(G,1,&ind,&val,ADD_VALUES);
359: }
360: }
361: }
363: /* Compute gradient contributions over the upper triangular elements */
364: for (j=0; j<=ny; j++) {
365: for (i=0; i<=nx; i++) {
366: k = nx*j + i;
367: vb = zero;
368: vl = zero;
369: v = zero;
370: if (i < nx && j > 0) vb = x[k-nx];
371: if (i > 0 && j < ny) vl = x[k-1];
372: if (i < nx && j < ny) v = x[k];
373: dvdx = (v-vl)/hx;
374: dvdy = (v-vb)/hy;
375: if (i != nx && j != 0) {
376: ind = k-nx; val = - dvdy/hy - cdiv3;
377: VecSetValues(G,1,&ind,&val,ADD_VALUES);
378: }
379: if (i != 0 && j != ny) {
380: ind = k-1; val = - dvdx/hx - cdiv3;
381: VecSetValues(G,1,&ind,&val,ADD_VALUES);
382: }
383: if (i != nx && j != ny) {
384: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
385: VecSetValues(G,1,&ind,&val,ADD_VALUES);
386: }
387: }
388: }
389: VecRestoreArrayRead(X,&x);
391: /* Assemble gradient vector */
392: VecAssemblyBegin(G);
393: VecAssemblyEnd(G);
395: /* Scale the gradient */
396: area = p5*hx*hy;
397: VecScale(G, area);
398: PetscLogFlops(24.0*nx*ny);
399: return 0;
400: }
402: /* ------------------------------------------------------------------- */
403: /*
404: FormHessian - Forms the Hessian matrix.
406: Input Parameters:
407: . tao - the Tao context
408: . X - the input vector
409: . ptr - optional user-defined context, as set by TaoSetHessian()
411: Output Parameters:
412: . H - Hessian matrix
413: . PrecH - optionally different preconditioning Hessian
414: . flag - flag indicating matrix structure
416: Notes:
417: This routine is intended simply as an example of the interface
418: to a Hessian evaluation routine. Since this example compute the
419: Hessian a column at a time, it is not particularly efficient and
420: is not recommended.
421: */
422: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
423: {
424: AppCtx *user = (AppCtx *) ptr;
425: PetscInt i,j, ndim = user->ndim;
426: PetscReal *y, zero = 0.0, one = 1.0;
427: PetscBool assembled;
430: user->xvec = X;
432: /* Initialize Hessian entries and work vector to zero */
433: MatAssembled(H,&assembled);
434: if (assembled)MatZeroEntries(H);
436: VecSet(user->s, zero);
438: /* Loop over matrix columns to compute entries of the Hessian */
439: for (j=0; j<ndim; j++) {
440: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
441: VecAssemblyBegin(user->s);
442: VecAssemblyEnd(user->s);
444: HessianProduct(ptr,user->s,user->y);
446: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
447: VecAssemblyBegin(user->s);
448: VecAssemblyEnd(user->s);
450: VecGetArray(user->y,&y);
451: for (i=0; i<ndim; i++) {
452: if (y[i] != zero) {
453: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
454: }
455: }
456: VecRestoreArray(user->y,&y);
457: }
458: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
459: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
460: return 0;
461: }
463: /* ------------------------------------------------------------------- */
464: /*
465: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
466: products.
468: Input Parameters:
469: . tao - the Tao context
470: . X - the input vector
471: . ptr - optional user-defined context, as set by TaoSetHessian()
473: Output Parameters:
474: . H - Hessian matrix
475: . PrecH - optionally different preconditioning Hessian
476: . flag - flag indicating matrix structure
477: */
478: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
479: {
480: AppCtx *user = (AppCtx *) ptr;
482: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
484: user->xvec = X;
485: return 0;
486: }
488: /* ------------------------------------------------------------------- */
489: /*
490: HessianProductMat - Computes the matrix-vector product
491: y = mat*svec.
493: Input Parameters:
494: . mat - input matrix
495: . svec - input vector
497: Output Parameters:
498: . y - solution vector
499: */
500: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
501: {
502: void *ptr;
505: MatShellGetContext(mat,&ptr);
506: HessianProduct(ptr,svec,y);
507: return 0;
508: }
510: /* ------------------------------------------------------------------- */
511: /*
512: Hessian Product - Computes the matrix-vector product:
513: y = f''(x)*svec.
515: Input Parameters:
516: . ptr - pointer to the user-defined context
517: . svec - input vector
519: Output Parameters:
520: . y - product vector
521: */
522: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
523: {
524: AppCtx *user = (AppCtx *)ptr;
525: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
526: const PetscScalar *x, *s;
527: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
528: PetscInt nx, ny, i, j, k, ind;
531: nx = user->mx;
532: ny = user->my;
533: hx = user->hx;
534: hy = user->hy;
535: hxhx = one/(hx*hx);
536: hyhy = one/(hy*hy);
538: /* Get pointers to vector data */
539: VecGetArrayRead(user->xvec,&x);
540: VecGetArrayRead(svec,&s);
542: /* Initialize product vector to zero */
543: VecSet(y, zero);
545: /* Compute f''(x)*s over the lower triangular elements */
546: for (j=-1; j<ny; j++) {
547: for (i=-1; i<nx; i++) {
548: k = nx*j + i;
549: v = zero;
550: vr = zero;
551: vt = zero;
552: if (i != -1 && j != -1) v = s[k];
553: if (i != nx-1 && j != -1) {
554: vr = s[k+1];
555: ind = k+1; val = hxhx*(vr-v);
556: VecSetValues(y,1,&ind,&val,ADD_VALUES);
557: }
558: if (i != -1 && j != ny-1) {
559: vt = s[k+nx];
560: ind = k+nx; val = hyhy*(vt-v);
561: VecSetValues(y,1,&ind,&val,ADD_VALUES);
562: }
563: if (i != -1 && j != -1) {
564: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
565: VecSetValues(y,1,&ind,&val,ADD_VALUES);
566: }
567: }
568: }
570: /* Compute f''(x)*s over the upper triangular elements */
571: for (j=0; j<=ny; j++) {
572: for (i=0; i<=nx; i++) {
573: k = nx*j + i;
574: v = zero;
575: vl = zero;
576: vb = zero;
577: if (i != nx && j != ny) v = s[k];
578: if (i != nx && j != 0) {
579: vb = s[k-nx];
580: ind = k-nx; val = hyhy*(vb-v);
581: VecSetValues(y,1,&ind,&val,ADD_VALUES);
582: }
583: if (i != 0 && j != ny) {
584: vl = s[k-1];
585: ind = k-1; val = hxhx*(vl-v);
586: VecSetValues(y,1,&ind,&val,ADD_VALUES);
587: }
588: if (i != nx && j != ny) {
589: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
590: VecSetValues(y,1,&ind,&val,ADD_VALUES);
591: }
592: }
593: }
594: /* Restore vector data */
595: VecRestoreArrayRead(svec,&s);
596: VecRestoreArrayRead(user->xvec,&x);
598: /* Assemble vector */
599: VecAssemblyBegin(y);
600: VecAssemblyEnd(y);
602: /* Scale resulting vector by area */
603: area = p5*hx*hy;
604: VecScale(y, area);
605: PetscLogFlops(18.0*nx*ny);
606: return 0;
607: }
609: /*TEST
611: build:
612: requires: !complex
614: test:
615: suffix: 1
616: args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4
618: test:
619: suffix: 2
620: args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4
622: test:
623: suffix: 3
624: args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
626: test:
627: suffix: 4
628: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
630: test:
631: suffix: 5
632: args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
634: test:
635: suffix: 6
636: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
638: TEST*/