Actual source code: eptorsion2.c
1: /* Program usage: mpirun -np <proc> eptorsion2 [-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 uniprocessor version of this code is eptorsion1.c; the Fortran
16: version of this code is eptorsion2f.F.
18: This application solves the problem without calculating hessians
19: ---------------------------------------------------------------------- */
21: /*
22: Include "tao.h" so that we can use TAO solvers. Note that this
23: file automatically includes files for lower-level support, such as those
24: provided by the PETSc library:
25: petsc.h - base PETSc routines petscvec.h - vectors
26: petscsys.h - sysem routines petscmat.h - matrices
27: petscis.h - index sets petscksp.h - Krylov subspace methods
28: petscviewer.h - viewers petscpc.h - preconditioners
29: Include "petscdm.h" so that we can use distributed arrays (DMs) for managing
30: the parallel mesh.
31: */
33: #include "taosolver.h"
34: #include "petscdm.h"
36: static char help[] =
37: "Demonstrates use of the TAO package to solve \n\
38: unconstrained minimization problems in parallel. This example is based on \n\
39: the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.\n\
40: The command line options are:\n\
41: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
42: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
43: -par <param>, where <param> = angle of twist per unit length\n\n";
45: /*T
46: Concepts: TAO - Solving an unconstrained minimization problem
47: Routines: TaoInitialize(); TaoFinalize();
48: Routines: TaoCreate(); TaoSetType();
49: Routines: TaoSetInitialVector();
50: Routines: TaoSetObjectiveAndGradientRoutine();
51: Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
52: Routines: TaoSolve();
53: Routines: TaoGetTerminationReason(); TaoDestroy();
54: Processors: n
55: T*/
57: /*
58: User-defined application context - contains data needed by the
59: application-provided call-back routines, FormFunction() and
60: FormGradient().
61: */
62: typedef struct {
63: /* parameters */
64: PetscInt mx, my; /* global discretization in x- and y-directions */
65: PetscReal param; /* nonlinearity parameter */
67: /* work space */
68: Vec localX; /* local vectors */
69: DM dm; /* distributed array data structure */
70: } AppCtx;
73: PetscErrorCode FormInitialGuess(AppCtx*, Vec);
74: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal*,Vec,void*);
75: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);
80: int main(int argc, char **argv)
81: {
83: Vec x;
84: Mat H;
85: PetscInt Nx, Ny;
86: TaoSolver tao;
87: TaoSolverTerminationReason reason;
88: PetscBool flg;
89: AppCtx user;
91: /* Initialize PETSc, TAO */
92: PetscInitialize(&argc, &argv, (char *)0, help);
93: TaoInitialize(&argc, &argv, (char *)0, help);
95: /* Specify default dimension of the problem */
96: user.param = 5.0; user.mx = 10; user.my = 10;
97: Nx = Ny = PETSC_DECIDE;
99: /* Check for any command line arguments that override defaults */
100: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,&flg);
101: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg);
102: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg);
104: PetscPrintf(PETSC_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n");
105: PetscPrintf(PETSC_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
107: /* Set up distributed array */
108: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,
109: DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
110: user.mx,user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
111: &user.dm);
113: /* Create vectors */
114: DMCreateGlobalVector(user.dm,&x);
116: DMCreateLocalVector(user.dm,&user.localX);
118: /* Create Hessian */
119: DMGetMatrix(user.dm,MATAIJ,&H);
120: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
122: /* The TAO code begins here */
124: /* Create TAO solver and set desired solution method */
125: TaoCreate(PETSC_COMM_WORLD,&tao);
126: TaoSetType(tao,"tao_cg");
128: /* Set initial solution guess */
129: FormInitialGuess(&user,x);
130: TaoSetInitialVector(tao,x);
132: /* Set routine for function and gradient evaluation */
133: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
135: TaoSetHessianRoutine(tao,H,H,FormHessian,(void*)&user);
138: /* Check for any TAO command line options */
139: TaoSetFromOptions(tao);
141: /* SOLVE THE APPLICATION */
142: TaoSolve(tao);
144: /* Get information on termination */
145: TaoGetTerminationReason(tao,&reason);
146: if (reason <= 0){
147: ierr=PetscPrintf(MPI_COMM_WORLD, "Try another method! \n");
148:
149: }
151: /* Free TAO data structures */
152: TaoDestroy(&tao);
154: /* Free PETSc data structures */
155: VecDestroy(&x);
156: MatDestroy(&H);
158: VecDestroy(&user.localX);
159: DMDestroy(&user.dm);
162: /* Finalize TAO and PETSc*/
163: TaoFinalize();
164: PetscFinalize();
165:
166: return 0;
167: }
170: /* ------------------------------------------------------------------- */
173: /*
174: FormInitialGuess - Computes an initial approximation to the solution.
176: Input Parameters:
177: . user - user-defined application context
178: . X - vector
180: Output Parameters:
181: X - vector
182: */
183: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
184: {
185: PetscErrorCode ierr;
186: PetscInt i, j, k, mx = user->mx, my = user->my;
187: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
188: PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;
191: /* Get local mesh boundaries */
192: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
193: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
195: /* Compute initial guess over locally owned part of mesh */
196: xe = xs+xm;
197: ye = ys+ym;
198: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
199: temp = PetscMin(j+1,my-j)*hy;
200: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
201: k = (j-gys)*gxm + i-gxs;
202: val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
203: VecSetValuesLocal(X,1,&k,&val,ADD_VALUES);
204: }
205: }
206: VecAssemblyBegin(X);
207: VecAssemblyEnd(X);
208: return(0);
209: }
212: /* ------------------------------------------------------------------ */
215: /*
216: FormFunctionGradient - Evaluates the function and corresponding gradient.
217:
218: Input Parameters:
219: tao - the TaoSolver context
220: X - the input vector
221: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
223: Output Parameters:
224: f - the newly evaluated function
225: G - the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(TaoSolver tao,Vec X,PetscReal *f,Vec G,void *ptr){
229: AppCtx *user = (AppCtx *)ptr;
230: PetscErrorCode ierr;
231: PetscInt i,j,k,ind;
232: PetscInt xe,ye,xsm,ysm,xep,yep;
233: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
234: PetscInt mx = user->mx, my = user->my;
235: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
236: PetscReal p5 = 0.5, area, val, flin, fquad;
237: PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
238: PetscReal hx = 1.0/(user->mx + 1);
239: PetscReal hy = 1.0/(user->my + 1);
240: Vec localX = user->localX;
244: /* Initialize */
245: flin = fquad = zero;
247: VecSet(G, zero);
248: /*
249: Scatter ghost points to local vector,using the 2-step process
250: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
251: By placing code between these two statements, computations can be
252: done while messages are in transition.
253: */
254: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
255: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
257: /* Get pointer to vector data */
258: VecGetArray(localX,&x);
260: /* Get local mesh boundaries */
261: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
262: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
264: /* Set local loop dimensions */
265: xe = xs+xm;
266: ye = ys+ym;
267: if (xs == 0) xsm = xs-1;
268: else xsm = xs;
269: if (ys == 0) ysm = ys-1;
270: else ysm = ys;
271: if (xe == mx) xep = xe+1;
272: else xep = xe;
273: if (ye == my) yep = ye+1;
274: else yep = ye;
276: /* Compute local gradient contributions over the lower triangular elements */
277: for (j=ysm; j<ye; j++) { /* for (j=-1; j<my; j++) */
278: for (i=xsm; i<xe; i++) { /* for (i=-1; i<mx; i++) */
279: k = (j-gys)*gxm + i-gxs;
280: v = zero;
281: vr = zero;
282: vt = zero;
283: if (i >= 0 && j >= 0) v = x[k];
284: if (i < mx-1 && j > -1) vr = x[k+1];
285: if (i > -1 && j < my-1) vt = x[k+gxm];
286: dvdx = (vr-v)/hx;
287: dvdy = (vt-v)/hy;
288: if (i != -1 && j != -1) {
289: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
290: VecSetValuesLocal(G,1,&k,&val,ADD_VALUES);
291: }
292: if (i != mx-1 && j != -1) {
293: ind = k+1; val = dvdx/hx - cdiv3;
294: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
295: }
296: if (i != -1 && j != my-1) {
297: ind = k+gxm; val = dvdy/hy - cdiv3;
298: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
299: }
300: fquad += dvdx*dvdx + dvdy*dvdy;
301: flin -= cdiv3 * (v + vr + vt);
302: }
303: }
305: /* Compute local gradient contributions over the upper triangular elements */
306: for (j=ys; j<yep; j++) { /* for (j=0; j<=my; j++) */
307: for (i=xs; i<xep; i++) { /* for (i=0; i<=mx; i++) */
308: k = (j-gys)*gxm + i-gxs;
309: vb = zero;
310: vl = zero;
311: v = zero;
312: if (i < mx && j > 0) vb = x[k-gxm];
313: if (i > 0 && j < my) vl = x[k-1];
314: if (i < mx && j < my) v = x[k];
315: dvdx = (v-vl)/hx;
316: dvdy = (v-vb)/hy;
317: if (i != mx && j != 0) {
318: ind = k-gxm; val = - dvdy/hy - cdiv3;
319: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
320: }
321: if (i != 0 && j != my) {
322: ind = k-1; val = - dvdx/hx - cdiv3;
323: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
324: }
325: if (i != mx && j != my) {
326: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
327: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
328: }
329: fquad += dvdx*dvdx + dvdy*dvdy;
330: flin -= cdiv3 * (vb + vl + v);
331: }
332: }
335: /* Restore vector */
336: VecRestoreArray(localX,&x);
338: /* Assemble gradient vector */
339: VecAssemblyBegin(G);
340: VecAssemblyEnd(G);
342: /* Scale the gradient */
343: area = p5*hx*hy;
344: floc = area * (p5 * fquad + flin);
345: VecScale(G, area);
347: /* Sum function contributions from all processes */
348: (PetscErrorCode)MPI_Allreduce((void*)&floc,(void*)f,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
350: ierr=PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16);
352: return(0);
353: }
359: PetscErrorCode FormHessian(TaoSolver tao, Vec X, Mat *H, Mat *Hpre, MatStructure *flag, void*ctx){
361: AppCtx *user= (AppCtx*) ctx;
363: PetscInt i,j,k;
364: PetscInt col[5],row;
365: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
366: PetscReal v[5];
367: 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;
368: Mat A=*H;
370: /* Compute the quadratic term in the objective function */
372: /*
373: Get local grid boundaries
374: */
377: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
378: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
380: for (j=ys; j<ys+ym; j++){
381:
382: for (i=xs; i< xs+xm; i++){
384: row=(j-gys)*gxm + (i-gxs);
386: k=0;
387: if (j>gys){
388: v[k]=-2*hyhy; col[k]=row - gxm; k++;
389: }
391: if (i>gxs){
392: v[k]= -2*hxhx; col[k]=row - 1; k++;
393: }
395: v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;
397: if (i+1 < gxs+gxm){
398: v[k]= -2.0*hxhx; col[k]=row+1; k++;
399: }
401: if (j+1 <gys+gym){
402: v[k]= -2*hyhy; col[k] = row+gxm; k++;
403: }
405: MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES);
407: }
408: }
409: /*
410: Assemble matrix, using the 2-step process:
411: MatAssemblyBegin(), MatAssemblyEnd().
412: By placing code between these two statements, computations can be
413: done while messages are in transition.
414: */
415: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
416: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
417: /*
418: Tell the matrix we will never add a new nonzero location to the
419: matrix. If we do it will generate an error.
420: */
421: MatScale(A,area);
422: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
423: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
425: PetscLogFlops(9*xm*ym+49*xm);
427: return(0);
428: }