Actual source code: eptorsion2.c
1: /*$Id$*/
3: /* Program usage: mpirun -np <proc> eptorsion2 [-help] [all TAO options] */
5: /* ----------------------------------------------------------------------
7: Elastic-plastic torsion problem.
9: The elastic plastic torsion problem arises from the determination
10: of the stress field on an infinitely long cylindrical bar, which is
11: equivalent to the solution of the following problem:
13: min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
14:
15: where C is the torsion angle per unit length.
17: The uniprocessor version of this code is eptorsion1.c; the Fortran
18: version of this code is eptorsion2f.F.
20: This application solves the problem without calculating hessians
21: ---------------------------------------------------------------------- */
23: /*
24: Include "tao.h" so that we can use TAO solvers. Note that this
25: file automatically includes files for lower-level support, such as those
26: provided by the PETSc library:
27: petsc.h - base PETSc routines petscvec.h - vectors
28: petscsys.h - sysem routines petscmat.h - matrices
29: petscis.h - index sets petscksp.h - Krylov subspace methods
30: petscviewer.h - viewers petscpc.h - preconditioners
31: Include "petscda.h" so that we can use distributed arrays (DAs) for managing
32: the parallel mesh.
33: */
35: #include "tao.h"
36: #include "petscda.h"
38: static char help[] =
39: "Demonstrates use of the TAO package to solve \n\
40: unconstrained minimization problems in parallel. This example is based on \n\
41: the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.\n\
42: The command line options are:\n\
43: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
44: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
45: -par <param>, where <param> = angle of twist per unit length\n\n";
47: /*T
48: Concepts: TAO - Solving an unconstrained minimization problem
49: Routines: TaoInitialize(); TaoFinalize();
50: Routines: TaoApplicationCreate(); TaoAppDestroy();
51: Routines: TaoCreate(); TaoDestroy();
52: Routines: TaoAppSetObjectiveAndGradientRoutine();
53: Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
54: Routines: TaoSetOptions();
55: Routines: TaoAppSetInitialSolutionVec();
56: Routines: TaoSolve(); TaoView(); TaoAppGetKSP();
57: Routines: TaoGetSolutionStatus();
58: Processors: n
59: T*/
63: /*
64: User-defined application context - contains data needed by the
65: application-provided call-back routines, FormFunction() and
66: FormGradient().
67: */
68: typedef struct {
69: /* parameters */
70: PetscInt mx, my; /* global discretization in x- and y-directions */
71: PetscReal param; /* nonlinearity parameter */
73: /* work space */
74: Vec localX; /* local vectors */
75: DA da; /* distributed array data structure */
76: } AppCtx;
78: /* -------- User-defined Routines --------- */
80: int FormInitialGuess(AppCtx*,Vec);
81: int FormFunctionGradient(TAO_APPLICATION,Vec,double*,Vec,void*);
82: int ComputeHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*);
86: int main(int argc,char **argv)
87: {
88: int info; /* used to check for functions returning nonzeros */
89: Vec x; /* solution, gradient vectors */
90: Mat H; /* Hessian matrix */
91: PetscInt Nx, Ny; /* number of processes in x- and y-directions */
92: TAO_SOLVER tao; /* TAO_SOLVER solver context */
93: TAO_APPLICATION torsionapp; /* TAO application context */
94: TaoTerminateReason reason;
95: PetscTruth flg;
96: PetscInt iter; /* iteration information */
97: double ff,gnorm;
98: AppCtx user; /* application context */
99: KSP ksp;
101: /* Initialize TAO, PETSc contexts */
102: info = PetscInitialize(&argc,&argv,(char *)0,help);
103: info = TaoInitialize(&argc,&argv,(char *)0,help);
105: /* Specify default parameters */
106: user.param = 5.0; user.mx = 10; user.my = 10;
107: Nx = Ny = PETSC_DECIDE;
109: /* Check for any command line arguments that override defaults */
110: info = PetscOptionsGetReal(TAO_NULL,"-par",&user.param,&flg); CHKERRQ(info);
111: info = PetscOptionsGetInt(TAO_NULL,"-my",&user.my,&flg); CHKERRQ(info);
112: info = PetscOptionsGetInt(TAO_NULL,"-mx",&user.mx,&flg); CHKERRQ(info);
114: PetscPrintf(PETSC_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n");
115: PetscPrintf(PETSC_COMM_WORLD,"mx: %d my: %d \n\n",user.mx,user.my);
117: /* Set up distributed array */
118: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,
119: user.mx,user.my,Nx,Ny,1,1,TAO_NULL,TAO_NULL,
120: &user.da); CHKERRQ(info);
122: /* Create vectors */
123: info = DACreateGlobalVector(user.da,&x); CHKERRQ(info);
125: info = DACreateLocalVector(user.da,&user.localX); CHKERRQ(info);
127: /* Create Hessian */
128: info = DAGetMatrix(user.da,MATAIJ,&H); CHKERRQ(info);
129: info = MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);
131: /* The TAO code begins here */
133: /* Create TAO solver and set desired solution method */
134: info = TaoCreate(MPI_COMM_WORLD,"tao_cg",&tao); CHKERRQ(info);
135: info = TaoApplicationCreate(PETSC_COMM_WORLD,&torsionapp); CHKERRQ(info);
137: /* Set initial solution guess */
138: info = FormInitialGuess(&user,x); CHKERRQ(info);
139: info = TaoAppSetInitialSolutionVec(torsionapp,x); CHKERRQ(info);
141: /* Set routine for function and gradient evaluation */
142: info = TaoAppSetObjectiveAndGradientRoutine(torsionapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);
144: info = TaoAppSetHessianMat(torsionapp,H,H); CHKERRQ(info);
145: info = TaoAppSetHessianRoutine(torsionapp,ComputeHessian,(void*)&user); CHKERRQ(info);
147: info = TaoAppGetKSP(torsionapp,&ksp); CHKERRQ(info);
148: if (ksp) { /* Modify the PETSc KSP structure */
149: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
150: }
152: /* Check for any TAO command line options */
153: info = TaoSetOptions(torsionapp,tao); CHKERRQ(info);
155: /* SOLVE THE APPLICATION */
156: info = TaoSolveApplication(torsionapp,tao); CHKERRQ(info);
158: /* Get information on termination */
159: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
160: if (reason <= 0){
161: PetscPrintf(MPI_COMM_WORLD, "Try another method! Iterations: %d, f: %4.2e, residual: %4.2e\n",
162: iter,ff,gnorm); CHKERRQ(info);
163: }
165: /*
166: To View TAO solver information use
167: info = TaoView(tao); CHKERRQ(info);
168: */
170: /* Free TAO data structures */
171: info = TaoDestroy(tao); CHKERRQ(info);
172: info = TaoAppDestroy(torsionapp); CHKERRQ(info);
174: /* Free PETSc data structures */
175: info = VecDestroy(x); CHKERRQ(info);
176: info = MatDestroy(H); CHKERRQ(info);
178: info = VecDestroy(user.localX); CHKERRQ(info);
179: info = DADestroy(user.da); CHKERRQ(info);
182: /* Finalize TAO and PETSc*/
183: TaoFinalize();
184: PetscFinalize();
186: return 0;
187: }
190: /* ------------------------------------------------------------------- */
193: /*
194: FormInitialGuess - Computes an initial approximation to the solution.
196: Input Parameters:
197: . user - user-defined application context
198: . X - vector
200: Output Parameters:
201: X - vector
202: */
203: int FormInitialGuess(AppCtx *user,Vec X)
204: {
205: int info;
206: PetscInt i, j, k, mx = user->mx, my = user->my;
207: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
208: PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;
210: /* Get local mesh boundaries */
211: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
212: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
214: /* Compute initial guess over locally owned part of mesh */
215: xe = xs+xm;
216: ye = ys+ym;
217: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
218: temp = TaoMin(j+1,my-j)*hy;
219: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
220: k = (j-gys)*gxm + i-gxs;
221: val = TaoMin((TaoMin(i+1,mx-i))*hx,temp);
222: info = VecSetValuesLocal(X,1,&k,&val,ADD_VALUES); CHKERRQ(info);
223: }
224: }
225: info = VecAssemblyBegin(X); CHKERRQ(info);
226: info = VecAssemblyEnd(X); CHKERRQ(info);
228: return 0;
229: }
232: /* ------------------------------------------------------------------ */
235: /*
236: FormFunctionGradient - Evaluates the function and corresponding gradient.
237:
238: Input Parameters:
239: taoapp - the TAO_APPLICATION context
240: X - the input vector
241: ptr - optional user-defined context, as set by TaoSetFunction()
243: Output Parameters:
244: f - the newly evaluated function
245: G - the newly evaluated gradient
246: */
247: int FormFunctionGradient(TAO_APPLICATION taoapp,Vec X,double *f,Vec G,void *ptr){
249: AppCtx *user = (AppCtx *)ptr;
250: int info;
251: PetscInt i,j,k,ind;
252: PetscInt xe,ye,xsm,ysm,xep,yep;
253: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
254: PetscInt mx = user->mx, my = user->my;
255: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
256: PetscReal p5 = 0.5, area, val, flin, fquad;
257: PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
258: PetscReal hx = 1.0/(user->mx + 1);
259: PetscReal hy = 1.0/(user->my + 1);
260: Vec localX = user->localX;
263: /* Initialize */
264: flin = fquad = zero;
266: info = VecSet(G, zero); CHKERRQ(info);
267: /*
268: Scatter ghost points to local vector,using the 2-step process
269: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
270: By placing code between these two statements, computations can be
271: done while messages are in transition.
272: */
273: info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
274: info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
276: /* Get pointer to vector data */
277: info = VecGetArray(localX,&x); CHKERRQ(info);
279: /* Get local mesh boundaries */
280: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
281: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
283: /* Set local loop dimensions */
284: xe = xs+xm;
285: ye = ys+ym;
286: if (xs == 0) xsm = xs-1;
287: else xsm = xs;
288: if (ys == 0) ysm = ys-1;
289: else ysm = ys;
290: if (xe == mx) xep = xe+1;
291: else xep = xe;
292: if (ye == my) yep = ye+1;
293: else yep = ye;
295: /* Compute local gradient contributions over the lower triangular elements */
296: for (j=ysm; j<ye; j++) { /* for (j=-1; j<my; j++) */
297: for (i=xsm; i<xe; i++) { /* for (i=-1; i<mx; i++) */
298: k = (j-gys)*gxm + i-gxs;
299: v = zero;
300: vr = zero;
301: vt = zero;
302: if (i >= 0 && j >= 0) v = x[k];
303: if (i < mx-1 && j > -1) vr = x[k+1];
304: if (i > -1 && j < my-1) vt = x[k+gxm];
305: dvdx = (vr-v)/hx;
306: dvdy = (vt-v)/hy;
307: if (i != -1 && j != -1) {
308: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
309: info = VecSetValuesLocal(G,1,&k,&val,ADD_VALUES); CHKERRQ(info);
310: }
311: if (i != mx-1 && j != -1) {
312: ind = k+1; val = dvdx/hx - cdiv3;
313: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
314: }
315: if (i != -1 && j != my-1) {
316: ind = k+gxm; val = dvdy/hy - cdiv3;
317: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
318: }
319: fquad += dvdx*dvdx + dvdy*dvdy;
320: flin -= cdiv3 * (v + vr + vt);
321: }
322: }
324: /* Compute local gradient contributions over the upper triangular elements */
325: for (j=ys; j<yep; j++) { /* for (j=0; j<=my; j++) */
326: for (i=xs; i<xep; i++) { /* for (i=0; i<=mx; i++) */
327: k = (j-gys)*gxm + i-gxs;
328: vb = zero;
329: vl = zero;
330: v = zero;
331: if (i < mx && j > 0) vb = x[k-gxm];
332: if (i > 0 && j < my) vl = x[k-1];
333: if (i < mx && j < my) v = x[k];
334: dvdx = (v-vl)/hx;
335: dvdy = (v-vb)/hy;
336: if (i != mx && j != 0) {
337: ind = k-gxm; val = - dvdy/hy - cdiv3;
338: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
339: }
340: if (i != 0 && j != my) {
341: ind = k-1; val = - dvdx/hx - cdiv3;
342: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
343: }
344: if (i != mx && j != my) {
345: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
346: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
347: }
348: fquad += dvdx*dvdx + dvdy*dvdy;
349: flin -= cdiv3 * (vb + vl + v);
350: }
351: }
354: /* Restore vector */
355: info = VecRestoreArray(localX,&x); CHKERRQ(info);
357: /* Assemble gradient vector */
358: info = VecAssemblyBegin(G); CHKERRQ(info);
359: info = VecAssemblyEnd(G); CHKERRQ(info);
361: /* Scale the gradient */
362: area = p5*hx*hy;
363: floc = area * (p5 * fquad + flin);
364: info = VecScale(G, area); CHKERRQ(info);
366: /* Sum function contributions from all processes */
367: MPI_Allreduce((void*)&floc,(void*)f,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
369: info=PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16); CHKERRQ(info);
370:
371: return 0;
372: }
378: int ComputeHessian(TAO_APPLICATION taoapp, Vec X, Mat *H, Mat *Hpre, MatStructure *flag, void*ctx){
380: AppCtx *user= (AppCtx*) ctx;
381: int info;
382: PetscInt i,j,k;
383: PetscInt col[5],row;
384: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
385: PetscReal v[5];
386: PetscReal hx=1.0/(user->mx+1), hy=1.0/(user->my+1), hxhx=1.0/(hx*hx), hyhy=1.0/(hy*hy), area=0.5*hx*hy;
387: Mat A=*H;
389: /* Compute the quadratic term in the objective function */
391: /*
392: Get local grid boundaries
393: */
395: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
396: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
398: for (j=ys; j<ys+ym; j++){
399:
400: for (i=xs; i< xs+xm; i++){
402: row=(j-gys)*gxm + (i-gxs);
404: k=0;
405: if (j>gys){
406: v[k]=-2*hyhy; col[k]=row - gxm; k++;
407: }
409: if (i>gxs){
410: v[k]= -2*hxhx; col[k]=row - 1; k++;
411: }
413: v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;
415: if (i+1 < gxs+gxm){
416: v[k]= -2.0*hxhx; col[k]=row+1; k++;
417: }
419: if (j+1 <gys+gym){
420: v[k]= -2*hyhy; col[k] = row+gxm; k++;
421: }
423: info = MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);
425: }
426: }
427: /*
428: Assemble matrix, using the 2-step process:
429: MatAssemblyBegin(), MatAssemblyEnd().
430: By placing code between these two statements, computations can be
431: done while messages are in transition.
432: */
433: info = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
434: info = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
435: /*
436: Tell the matrix we will never add a new nonzero location to the
437: matrix. If we do it will generate an error.
438: */
439: info = MatScale(A,area); CHKERRQ(info);
440: info = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); CHKERRQ(info);
441: info = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);
443: info = PetscLogFlops(9*xm*ym+49*xm); CHKERRQ(info);
445: return 0;
446: }