Actual source code: eptorsion1.c
petsc-3.13.6 2020-09-29
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 Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
55: Section 1.5 Writing Application Codes with PETSc-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; /* Section 1.5 Writing Application Codes with PETSc context */
88: PetscMPIInt size; /* number of processes */
89: PetscReal one=1.0;
90:
91: PetscBool test_lmvm = PETSC_FALSE;
92: KSP ksp;
93: PC pc;
94: Mat M;
95: Vec in, out, out2;
96: PetscReal mult_solve_dist;
98: /* Initialize TAO,PETSc */
99: PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
100: MPI_Comm_size(MPI_COMM_WORLD,&size);
101: if (size >1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors");
103: /* Specify default parameters for the problem, check for command-line overrides */
104: user.param = 5.0;
105: PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
106: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
107: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
108: PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);
110: PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
111: PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my);
112: user.ndim = mx * my; user.mx = mx; user.my = my;
113: user.hx = one/(mx+1); user.hy = one/(my+1);
115: /* Allocate vectors */
116: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
117: VecDuplicate(user.y,&user.s);
118: VecDuplicate(user.y,&x);
120: /* Create TAO solver and set desired solution method */
121: TaoCreate(PETSC_COMM_SELF,&tao);
122: TaoSetType(tao,TAOLMVM);
124: /* Set solution vector with an initial guess */
125: FormInitialGuess(&user,x);
126: TaoSetInitialVector(tao,x);
128: /* Set routine for function and gradient evaluation */
129: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
131: /* From command line options, determine if using matrix-free hessian */
132: PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
133: if (flg) {
134: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
135: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
136: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
138: TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
139: } else {
140: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
141: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
142: TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
143: }
145: /* Test the LMVM matrix */
146: if (test_lmvm) {
147: PetscOptionsSetValue(NULL, "-tao_type", "bntr");
148: PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
149: }
151: /* Check for any TAO command line options */
152: TaoSetFromOptions(tao);
154: /* SOLVE THE APPLICATION */
155: TaoSolve(tao);
156:
157: /* Test the LMVM matrix */
158: if (test_lmvm) {
159: TaoGetKSP(tao, &ksp);
160: KSPGetPC(ksp, &pc);
161: PCLMVMGetMatLMVM(pc, &M);
162: VecDuplicate(x, &in);
163: VecDuplicate(x, &out);
164: VecDuplicate(x, &out2);
165: VecSet(in, 5.0);
166: MatMult(M, in, out);
167: MatSolve(M, out, out2);
168: VecAXPY(out2, -1.0, in);
169: VecNorm(out2, NORM_2, &mult_solve_dist);
170: PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist);
171: VecDestroy(&in);
172: VecDestroy(&out);
173: VecDestroy(&out2);
174: }
176: TaoDestroy(&tao);
177: VecDestroy(&user.s);
178: VecDestroy(&user.y);
179: VecDestroy(&x);
180: MatDestroy(&H);
182: PetscFinalize();
183: return ierr;
184: }
186: /* ------------------------------------------------------------------- */
187: /*
188: FormInitialGuess - Computes an initial approximation to the solution.
190: Input Parameters:
191: . user - user-defined Section 1.5 Writing Application Codes with PETSc context
192: . X - vector
194: Output Parameters:
195: . X - vector
196: */
197: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
198: {
199: PetscReal hx = user->hx, hy = user->hy, temp;
200: PetscReal val;
202: PetscInt i, j, k, nx = user->mx, ny = user->my;
204: /* Compute initial guess */
206: for (j=0; j<ny; j++) {
207: temp = PetscMin(j+1,ny-j)*hy;
208: for (i=0; i<nx; i++) {
209: k = nx*j + i;
210: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
211: VecSetValues(X,1,&k,&val,ADD_VALUES);
212: }
213: }
214: VecAssemblyBegin(X);
215: VecAssemblyEnd(X);
216: return(0);
217: }
219: /* ------------------------------------------------------------------- */
220: /*
221: FormFunctionGradient - Evaluates the function and corresponding gradient.
223: Input Parameters:
224: tao - the Tao context
225: X - the input vector
226: ptr - optional user-defined context, as set by TaoSetFunction()
228: Output Parameters:
229: f - the newly evaluated function
230: G - the newly evaluated gradient
231: */
232: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
233: {
237: FormFunction(tao,X,f,ptr);
238: FormGradient(tao,X,G,ptr);
239: return(0);
240: }
242: /* ------------------------------------------------------------------- */
243: /*
244: FormFunction - Evaluates the function, f(X).
246: Input Parameters:
247: . tao - the Tao context
248: . X - the input vector
249: . ptr - optional user-defined context, as set by TaoSetFunction()
251: Output Parameters:
252: . f - the newly evaluated function
253: */
254: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
255: {
256: AppCtx *user = (AppCtx *) ptr;
257: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
258: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
259: PetscReal v, cdiv3 = user->param/three;
260: const PetscScalar *x;
261: PetscErrorCode ierr;
262: PetscInt nx = user->mx, ny = user->my, i, j, k;
265: /* Get pointer to vector data */
266: VecGetArrayRead(X,&x);
268: /* Compute function contributions over the lower triangular elements */
269: for (j=-1; j<ny; j++) {
270: for (i=-1; i<nx; i++) {
271: k = nx*j + i;
272: v = zero;
273: vr = zero;
274: vt = zero;
275: if (i >= 0 && j >= 0) v = x[k];
276: if (i < nx-1 && j > -1) vr = x[k+1];
277: if (i > -1 && j < ny-1) vt = x[k+nx];
278: dvdx = (vr-v)/hx;
279: dvdy = (vt-v)/hy;
280: fquad += dvdx*dvdx + dvdy*dvdy;
281: flin -= cdiv3*(v+vr+vt);
282: }
283: }
285: /* Compute function contributions over the upper triangular elements */
286: for (j=0; j<=ny; j++) {
287: for (i=0; i<=nx; i++) {
288: k = nx*j + i;
289: vb = zero;
290: vl = zero;
291: v = zero;
292: if (i < nx && j > 0) vb = x[k-nx];
293: if (i > 0 && j < ny) vl = x[k-1];
294: if (i < nx && j < ny) v = x[k];
295: dvdx = (v-vl)/hx;
296: dvdy = (v-vb)/hy;
297: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
298: flin = flin - cdiv3*(vb+vl+v);
299: }
300: }
302: /* Restore vector */
303: VecRestoreArrayRead(X,&x);
305: /* Scale the function */
306: area = p5*hx*hy;
307: *f = area*(p5*fquad+flin);
309: PetscLogFlops(24.0*nx*ny);
310: return(0);
311: }
313: /* ------------------------------------------------------------------- */
314: /*
315: FormGradient - Evaluates the gradient, G(X).
317: Input Parameters:
318: . tao - the Tao context
319: . X - input vector
320: . ptr - optional user-defined context
322: Output Parameters:
323: . G - vector containing the newly evaluated gradient
324: */
325: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
326: {
327: AppCtx *user = (AppCtx *) ptr;
328: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val;
329: PetscErrorCode ierr;
330: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
331: PetscReal hx = user->hx, hy = user->hy;
332: PetscReal vb, vl, vr, vt, dvdx, dvdy;
333: PetscReal v, cdiv3 = user->param/three;
334: const PetscScalar *x;
337: /* Initialize gradient to zero */
338: VecSet(G, zero);
340: /* Get pointer to vector data */
341: VecGetArrayRead(X,&x);
343: /* Compute gradient contributions over the lower triangular elements */
344: for (j=-1; j<ny; j++) {
345: for (i=-1; i<nx; i++) {
346: k = nx*j + i;
347: v = zero;
348: vr = zero;
349: vt = zero;
350: if (i >= 0 && j >= 0) v = x[k];
351: if (i < nx-1 && j > -1) vr = x[k+1];
352: if (i > -1 && j < ny-1) vt = x[k+nx];
353: dvdx = (vr-v)/hx;
354: dvdy = (vt-v)/hy;
355: if (i != -1 && j != -1) {
356: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
357: VecSetValues(G,1,&ind,&val,ADD_VALUES);
358: }
359: if (i != nx-1 && j != -1) {
360: ind = k+1; val = dvdx/hx - cdiv3;
361: VecSetValues(G,1,&ind,&val,ADD_VALUES);
362: }
363: if (i != -1 && j != ny-1) {
364: ind = k+nx; val = dvdy/hy - cdiv3;
365: VecSetValues(G,1,&ind,&val,ADD_VALUES);
366: }
367: }
368: }
370: /* Compute gradient contributions over the upper triangular elements */
371: for (j=0; j<=ny; j++) {
372: for (i=0; i<=nx; i++) {
373: k = nx*j + i;
374: vb = zero;
375: vl = zero;
376: v = zero;
377: if (i < nx && j > 0) vb = x[k-nx];
378: if (i > 0 && j < ny) vl = x[k-1];
379: if (i < nx && j < ny) v = x[k];
380: dvdx = (v-vl)/hx;
381: dvdy = (v-vb)/hy;
382: if (i != nx && j != 0) {
383: ind = k-nx; val = - dvdy/hy - cdiv3;
384: VecSetValues(G,1,&ind,&val,ADD_VALUES);
385: }
386: if (i != 0 && j != ny) {
387: ind = k-1; val = - dvdx/hx - cdiv3;
388: VecSetValues(G,1,&ind,&val,ADD_VALUES);
389: }
390: if (i != nx && j != ny) {
391: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
392: VecSetValues(G,1,&ind,&val,ADD_VALUES);
393: }
394: }
395: }
396: VecRestoreArrayRead(X,&x);
398: /* Assemble gradient vector */
399: VecAssemblyBegin(G);
400: VecAssemblyEnd(G);
402: /* Scale the gradient */
403: area = p5*hx*hy;
404: VecScale(G, area);
405: PetscLogFlops(24.0*nx*ny);
406: return(0);
407: }
409: /* ------------------------------------------------------------------- */
410: /*
411: FormHessian - Forms the Hessian matrix.
413: Input Parameters:
414: . tao - the Tao context
415: . X - the input vector
416: . ptr - optional user-defined context, as set by TaoSetHessian()
418: Output Parameters:
419: . H - Hessian matrix
420: . PrecH - optionally different preconditioning Hessian
421: . flag - flag indicating matrix structure
423: Notes:
424: This routine is intended simply as an example of the interface
425: to a Hessian evaluation routine. Since this example compute the
426: Hessian a column at a time, it is not particularly efficient and
427: is not recommended.
428: */
429: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
430: {
431: AppCtx *user = (AppCtx *) ptr;
433: PetscInt i,j, ndim = user->ndim;
434: PetscReal *y, zero = 0.0, one = 1.0;
435: PetscBool assembled;
438: user->xvec = X;
440: /* Initialize Hessian entries and work vector to zero */
441: MatAssembled(H,&assembled);
442: if (assembled){MatZeroEntries(H); }
444: VecSet(user->s, zero);
446: /* Loop over matrix columns to compute entries of the Hessian */
447: for (j=0; j<ndim; j++) {
448: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
449: VecAssemblyBegin(user->s);
450: VecAssemblyEnd(user->s);
452: HessianProduct(ptr,user->s,user->y);
454: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
455: VecAssemblyBegin(user->s);
456: VecAssemblyEnd(user->s);
458: VecGetArray(user->y,&y);
459: for (i=0; i<ndim; i++) {
460: if (y[i] != zero) {
461: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
462: }
463: }
464: VecRestoreArray(user->y,&y);
465: }
466: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
467: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
468: return(0);
469: }
471: /* ------------------------------------------------------------------- */
472: /*
473: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
474: products.
476: Input Parameters:
477: . tao - the Tao context
478: . X - the input vector
479: . ptr - optional user-defined context, as set by TaoSetHessian()
481: Output Parameters:
482: . H - Hessian matrix
483: . PrecH - optionally different preconditioning Hessian
484: . flag - flag indicating matrix structure
485: */
486: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
487: {
488: AppCtx *user = (AppCtx *) ptr;
490: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
491: user->xvec = X;
492: return(0);
493: }
495: /* ------------------------------------------------------------------- */
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;
513: MatShellGetContext(mat,&ptr);
514: HessianProduct(ptr,svec,y);
515: return(0);
516: }
518: /* ------------------------------------------------------------------- */
519: /*
520: Hessian Product - Computes the matrix-vector product:
521: y = f''(x)*svec.
523: Input Parameters:
524: . ptr - pointer to the user-defined context
525: . svec - input vector
527: Output Parameters:
528: . y - product vector
529: */
530: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
531: {
532: AppCtx *user = (AppCtx *)ptr;
533: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
534: const PetscScalar *x, *s;
535: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
536: PetscErrorCode ierr;
537: PetscInt nx, ny, i, j, k, ind;
540: nx = user->mx;
541: ny = user->my;
542: hx = user->hx;
543: hy = user->hy;
544: hxhx = one/(hx*hx);
545: hyhy = one/(hy*hy);
547: /* Get pointers to vector data */
548: VecGetArrayRead(user->xvec,&x);
549: VecGetArrayRead(svec,&s);
551: /* Initialize product vector to zero */
552: VecSet(y, zero);
554: /* Compute f''(x)*s over the lower triangular elements */
555: for (j=-1; j<ny; j++) {
556: for (i=-1; i<nx; i++) {
557: k = nx*j + i;
558: v = zero;
559: vr = zero;
560: vt = zero;
561: if (i != -1 && j != -1) v = s[k];
562: if (i != nx-1 && j != -1) {
563: vr = s[k+1];
564: ind = k+1; val = hxhx*(vr-v);
565: VecSetValues(y,1,&ind,&val,ADD_VALUES);
566: }
567: if (i != -1 && j != ny-1) {
568: vt = s[k+nx];
569: ind = k+nx; val = hyhy*(vt-v);
570: VecSetValues(y,1,&ind,&val,ADD_VALUES);
571: }
572: if (i != -1 && j != -1) {
573: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
574: VecSetValues(y,1,&ind,&val,ADD_VALUES);
575: }
576: }
577: }
579: /* Compute f''(x)*s over the upper triangular elements */
580: for (j=0; j<=ny; j++) {
581: for (i=0; i<=nx; i++) {
582: k = nx*j + i;
583: v = zero;
584: vl = zero;
585: vb = zero;
586: if (i != nx && j != ny) v = s[k];
587: if (i != nx && j != 0) {
588: vb = s[k-nx];
589: ind = k-nx; val = hyhy*(vb-v);
590: VecSetValues(y,1,&ind,&val,ADD_VALUES);
591: }
592: if (i != 0 && j != ny) {
593: vl = s[k-1];
594: ind = k-1; val = hxhx*(vl-v);
595: VecSetValues(y,1,&ind,&val,ADD_VALUES);
596: }
597: if (i != nx && j != ny) {
598: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
599: VecSetValues(y,1,&ind,&val,ADD_VALUES);
600: }
601: }
602: }
603: /* Restore vector data */
604: VecRestoreArrayRead(svec,&s);
605: VecRestoreArrayRead(user->xvec,&x);
607: /* Assemble vector */
608: VecAssemblyBegin(y);
609: VecAssemblyEnd(y);
611: /* Scale resulting vector by area */
612: area = p5*hx*hy;
613: VecScale(y, area);
614: PetscLogFlops(18.0*nx*ny);
615: return(0);
616: }
619: /*TEST
621: build:
622: requires: !complex
624: test:
625: suffix: 1
626: args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4
628: test:
629: suffix: 2
630: args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4
632: test:
633: suffix: 3
634: args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
635:
636: test:
637: suffix: 4
638: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
639:
640: test:
641: suffix: 5
642: args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
643:
644: test:
645: suffix: 6
646: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
648: TEST*/