Actual source code: eptorsion2.c
1: /*$Id: eptorsion2.c 1.55 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/
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(); TaoGetKSP();
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: int 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: int 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: int 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: /* Set up distributed array */
115: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,
116: user.mx,user.my,Nx,Ny,1,1,TAO_NULL,TAO_NULL,
117: &user.da); CHKERRQ(info);
119: /* Create vectors */
120: info = DACreateGlobalVector(user.da,&x); CHKERRQ(info);
122: info = DACreateLocalVector(user.da,&user.localX); CHKERRQ(info);
124: /* Create Hessian */
125: info = DAGetMatrix(user.da,MATAIJ,&H); CHKERRQ(info);
126: info = MatSetOption(H,MAT_SYMMETRIC); CHKERRQ(info);
128: /* The TAO code begins here */
130: /* Create TAO solver and set desired solution method */
131: info = TaoCreate(MPI_COMM_WORLD,"tao_cg_fr",&tao); CHKERRQ(info);
132: info = TaoApplicationCreate(PETSC_COMM_WORLD,&torsionapp); CHKERRQ(info);
134: /* Set initial solution guess */
135: info = FormInitialGuess(&user,x); CHKERRQ(info);
136: info = TaoAppSetInitialSolutionVec(torsionapp,x); CHKERRQ(info);
138: /* Set routine for function and gradient evaluation */
139: info = TaoAppSetObjectiveAndGradientRoutine(torsionapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);
141: info = TaoAppSetHessianMat(torsionapp,H,H); CHKERRQ(info);
142: info = TaoAppSetHessianRoutine(torsionapp,ComputeHessian,(void*)&user); CHKERRQ(info);
144: info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
145: if (ksp) { /* Modify the PETSc KSP structure */
146: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
147: }
149: /* Check for any TAO command line options */
150: info = TaoSetOptions(torsionapp,tao); CHKERRQ(info);
152: /* SOLVE THE APPLICATION */
153: info = TaoSolveApplication(torsionapp,tao); CHKERRQ(info);
155: /* Get information on termination */
156: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
157: if (reason <= 0){
158: PetscPrintf(MPI_COMM_WORLD, "Try another method! Iterations: %d, f: %4.2e, residual: %4.2e\n",
159: iter,ff,gnorm); CHKERRQ(info);
160: }
162: /*
163: To View TAO solver information use
164: info = TaoView(tao); CHKERRQ(info);
165: */
167: /* Free TAO data structures */
168: info = TaoDestroy(tao); CHKERRQ(info);
169: info = TaoAppDestroy(torsionapp); CHKERRQ(info);
171: /* Free PETSc data structures */
172: info = VecDestroy(x); CHKERRQ(info);
173: info = MatDestroy(H); CHKERRQ(info);
175: info = VecDestroy(user.localX); CHKERRQ(info);
176: info = DADestroy(user.da); CHKERRQ(info);
179: /* Finalize TAO and PETSc*/
180: TaoFinalize();
181: PetscFinalize();
183: return 0;
184: }
187: /* ------------------------------------------------------------------- */
190: /*
191: FormInitialGuess - Computes an initial approximation to the solution.
193: Input Parameters:
194: . user - user-defined application context
195: . X - vector
197: Output Parameters:
198: X - vector
199: */
200: int FormInitialGuess(AppCtx *user,Vec X)
201: {
202: int info, i, j, k, mx = user->mx, my = user->my;
203: int xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
204: PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;
206: /* Get local mesh boundaries */
207: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
208: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
210: /* Compute initial guess over locally owned part of mesh */
211: xe = xs+xm;
212: ye = ys+ym;
213: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
214: temp = TaoMin(j+1,my-j)*hy;
215: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
216: k = (j-gys)*gxm + i-gxs;
217: val = TaoMin((TaoMin(i+1,mx-i))*hx,temp);
218: info = VecSetValuesLocal(X,1,&k,&val,ADD_VALUES); CHKERRQ(info);
219: }
220: }
222: return 0;
223: }
226: /* ------------------------------------------------------------------ */
229: /*
230: FormFunctionGradient - Evaluates the function and corresponding gradient.
231:
232: Input Parameters:
233: taoapp - the TAO_APPLICATION context
234: X - the input vector
235: ptr - optional user-defined context, as set by TaoSetFunction()
237: Output Parameters:
238: f - the newly evaluated function
239: G - the newly evaluated gradient
240: */
241: int FormFunctionGradient(TAO_APPLICATION taoapp,Vec X,double *f,Vec G,void *ptr){
243: AppCtx *user = (AppCtx *)ptr;
244: int info,i,j,k,ind;
245: int xe,ye,xsm,ysm,xep,yep;
246: int xs, ys, xm, ym, gxm, gym, gxs, gys;
247: int mx = user->mx, my = user->my;
248: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
249: PetscReal p5 = 0.5, area, val, flin, fquad;
250: PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
251: PetscReal hx = 1.0/(user->mx + 1);
252: PetscReal hy = 1.0/(user->my + 1);
253: Vec localX = user->localX;
256: /* Initialize */
257: flin = fquad = zero;
259: info = VecSet(G, zero); CHKERRQ(info);
260: /*
261: Scatter ghost points to local vector,using the 2-step process
262: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
263: By placing code between these two statements, computations can be
264: done while messages are in transition.
265: */
266: info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
267: info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
269: /* Get pointer to vector data */
270: info = VecGetArray(localX,&x); CHKERRQ(info);
272: /* Get local mesh boundaries */
273: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
274: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
276: /* Set local loop dimensions */
277: xe = xs+xm;
278: ye = ys+ym;
279: if (xs == 0) xsm = xs-1;
280: else xsm = xs;
281: if (ys == 0) ysm = ys-1;
282: else ysm = ys;
283: if (xe == mx) xep = xe+1;
284: else xep = xe;
285: if (ye == my) yep = ye+1;
286: else yep = ye;
288: /* Compute local gradient contributions over the lower triangular elements */
289: for (j=ysm; j<ye; j++) { /* for (j=-1; j<my; j++) */
290: for (i=xsm; i<xe; i++) { /* for (i=-1; i<mx; i++) */
291: k = (j-gys)*gxm + i-gxs;
292: v = zero;
293: vr = zero;
294: vt = zero;
295: if (i >= 0 && j >= 0) v = x[k];
296: if (i < mx-1 && j > -1) vr = x[k+1];
297: if (i > -1 && j < my-1) vt = x[k+gxm];
298: dvdx = (vr-v)/hx;
299: dvdy = (vt-v)/hy;
300: if (i != -1 && j != -1) {
301: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
302: info = VecSetValuesLocal(G,1,&k,&val,ADD_VALUES); CHKERRQ(info);
303: }
304: if (i != mx-1 && j != -1) {
305: ind = k+1; val = dvdx/hx - cdiv3;
306: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
307: }
308: if (i != -1 && j != my-1) {
309: ind = k+gxm; val = dvdy/hy - cdiv3;
310: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
311: }
312: fquad += dvdx*dvdx + dvdy*dvdy;
313: flin -= cdiv3 * (v + vr + vt);
314: }
315: }
317: /* Compute local gradient contributions over the upper triangular elements */
318: for (j=ys; j<yep; j++) { /* for (j=0; j<=my; j++) */
319: for (i=xs; i<xep; i++) { /* for (i=0; i<=mx; i++) */
320: k = (j-gys)*gxm + i-gxs;
321: vb = zero;
322: vl = zero;
323: v = zero;
324: if (i < mx && j > 0) vb = x[k-gxm];
325: if (i > 0 && j < my) vl = x[k-1];
326: if (i < mx && j < my) v = x[k];
327: dvdx = (v-vl)/hx;
328: dvdy = (v-vb)/hy;
329: if (i != mx && j != 0) {
330: ind = k-gxm; val = - dvdy/hy - cdiv3;
331: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
332: }
333: if (i != 0 && j != my) {
334: ind = k-1; val = - dvdx/hx - cdiv3;
335: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
336: }
337: if (i != mx && j != my) {
338: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
339: info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
340: }
341: fquad += dvdx*dvdx + dvdy*dvdy;
342: flin -= cdiv3 * (vb + vl + v);
343: }
344: }
347: /* Restore vector */
348: info = VecRestoreArray(localX,&x); CHKERRQ(info);
350: /* Assemble gradient vector */
351: info = VecAssemblyBegin(G); CHKERRQ(info);
352: info = VecAssemblyEnd(G); CHKERRQ(info);
354: /* Scale the gradient */
355: area = p5*hx*hy;
356: floc = area * (p5 * fquad + flin);
357: info = VecScale(G, area); CHKERRQ(info);
359: /* Sum function contributions from all processes */
360: MPI_Allreduce((void*)&floc,(void*)f,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
362: info=PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16); CHKERRQ(info);
363:
364: return 0;
365: }
371: int ComputeHessian(TAO_APPLICATION taoapp, Vec X, Mat *H, Mat *Hpre, MatStructure *flag, void*ctx){
373: AppCtx *user= (AppCtx*) ctx;
374: int i,j,k,info;
375: int col[5],row;
376: int xs,xm,gxs,gxm,ys,ym,gys,gym;
377: PetscReal v[5];
378: 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;
379: Mat A=*H;
381: /* Compute the quadratic term in the objective function */
383: /*
384: Get local grid boundaries
385: */
387: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
388: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
390: for (j=ys; j<ys+ym; j++){
391:
392: for (i=xs; i< xs+xm; i++){
394: row=(j-gys)*gxm + (i-gxs);
396: k=0;
397: if (j>gys){
398: v[k]=-2*hyhy; col[k]=row - gxm; k++;
399: }
401: if (i>gxs){
402: v[k]= -2*hxhx; col[k]=row - 1; k++;
403: }
405: v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;
407: if (i+1 < gxs+gxm){
408: v[k]= -2.0*hxhx; col[k]=row+1; k++;
409: }
411: if (j+1 <gys+gym){
412: v[k]= -2*hyhy; col[k] = row+gxm; k++;
413: }
415: info = MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);
417: }
418: }
419: /*
420: Assemble matrix, using the 2-step process:
421: MatAssemblyBegin(), MatAssemblyEnd().
422: By placing code between these two statements, computations can be
423: done while messages are in transition.
424: */
425: info = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
426: info = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
427: /*
428: Tell the matrix we will never add a new nonzero location to the
429: matrix. If we do it will generate an error.
430: */
431: info = MatScale(A,area); CHKERRQ(info);
432: info = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ(info);
433: info = MatSetOption(A,MAT_SYMMETRIC); CHKERRQ(info);
435: info = PetscLogFlops(9*xm*ym+49*xm); CHKERRQ(info);
437: return 0;
438: }