Actual source code: eptorsion2.c
1: /* Program usage: mpiexec -n <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)}
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 "petsctao.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 - system routines petscmat.h - matrices
27: petscis.h - index sets petscksp.h - Krylov subspace methods
28: petscviewer.h - viewers petscpc.h - preconditioners
29: Include "petscdmda.h" so that we can use distributed arrays (DMs) for managing
30: the parallel mesh.
31: */
33: #include <petsctao.h>
34: #include <petscdmda.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: /*
46: User-defined application context - contains data needed by the
47: application-provided call-back routines, FormFunction() and
48: FormGradient().
49: */
50: typedef struct {
51: /* parameters */
52: PetscInt mx, my; /* global discretization in x- and y-directions */
53: PetscReal param; /* nonlinearity parameter */
55: /* work space */
56: Vec localX; /* local vectors */
57: DM dm; /* distributed array data structure */
58: } AppCtx;
60: PetscErrorCode FormInitialGuess(AppCtx*, Vec);
61: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
62: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
64: int main(int argc, char **argv)
65: {
66: Vec x;
67: Mat H;
68: PetscInt Nx, Ny;
69: Tao tao;
70: PetscBool flg;
71: KSP ksp;
72: PC pc;
73: AppCtx user;
75: PetscInitialize(&argc, &argv, (char *)0, help);
77: /* Specify default dimension of the problem */
78: user.param = 5.0; user.mx = 10; user.my = 10;
79: Nx = Ny = PETSC_DECIDE;
81: /* Check for any command line arguments that override defaults */
82: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
83: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
84: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
86: PetscPrintf(PETSC_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n");
87: PetscPrintf(PETSC_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
89: /* Set up distributed array */
90: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
91: DMSetFromOptions(user.dm);
92: DMSetUp(user.dm);
94: /* Create vectors */
95: DMCreateGlobalVector(user.dm,&x);
97: DMCreateLocalVector(user.dm,&user.localX);
99: /* Create Hessian */
100: DMCreateMatrix(user.dm,&H);
101: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
103: /* The TAO code begins here */
105: /* Create TAO solver and set desired solution method */
106: TaoCreate(PETSC_COMM_WORLD,&tao);
107: TaoSetType(tao,TAOCG);
109: /* Set initial solution guess */
110: FormInitialGuess(&user,x);
111: TaoSetSolution(tao,x);
113: /* Set routine for function and gradient evaluation */
114: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
116: TaoSetHessian(tao,H,H,FormHessian,(void*)&user);
118: /* Check for any TAO command line options */
119: TaoSetFromOptions(tao);
121: TaoGetKSP(tao,&ksp);
122: if (ksp) {
123: KSPGetPC(ksp,&pc);
124: PCSetType(pc,PCNONE);
125: }
127: /* SOLVE THE APPLICATION */
128: TaoSolve(tao);
130: /* Free TAO data structures */
131: TaoDestroy(&tao);
133: /* Free PETSc data structures */
134: VecDestroy(&x);
135: MatDestroy(&H);
137: VecDestroy(&user.localX);
138: DMDestroy(&user.dm);
140: PetscFinalize();
141: return 0;
142: }
144: /* ------------------------------------------------------------------- */
145: /*
146: FormInitialGuess - Computes an initial approximation to the solution.
148: Input Parameters:
149: . user - user-defined application context
150: . X - vector
152: Output Parameters:
153: X - vector
154: */
155: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
156: {
157: PetscInt i, j, k, mx = user->mx, my = user->my;
158: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
159: PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;
161: /* Get local mesh boundaries */
162: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
163: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
165: /* Compute initial guess over locally owned part of mesh */
166: xe = xs+xm;
167: ye = ys+ym;
168: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
169: temp = PetscMin(j+1,my-j)*hy;
170: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
171: k = (j-gys)*gxm + i-gxs;
172: val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
173: VecSetValuesLocal(X,1,&k,&val,ADD_VALUES);
174: }
175: }
176: VecAssemblyBegin(X);
177: VecAssemblyEnd(X);
178: return 0;
179: }
181: /* ------------------------------------------------------------------ */
182: /*
183: FormFunctionGradient - Evaluates the function and corresponding gradient.
185: Input Parameters:
186: tao - the Tao context
187: X - the input vector
188: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
190: Output Parameters:
191: f - the newly evaluated function
192: G - the newly evaluated gradient
193: */
194: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
195: {
196: AppCtx *user = (AppCtx *)ptr;
197: PetscInt i,j,k,ind;
198: PetscInt xe,ye,xsm,ysm,xep,yep;
199: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
200: PetscInt mx = user->mx, my = user->my;
201: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
202: PetscReal p5 = 0.5, area, val, flin, fquad;
203: PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
204: PetscReal hx = 1.0/(user->mx + 1);
205: PetscReal hy = 1.0/(user->my + 1);
206: Vec localX = user->localX;
208: /* Initialize */
209: flin = fquad = zero;
211: VecSet(G, zero);
212: /*
213: Scatter ghost points to local vector,using the 2-step process
214: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
215: By placing code between these two statements, computations can be
216: done while messages are in transition.
217: */
218: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
219: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
221: /* Get pointer to vector data */
222: VecGetArray(localX,&x);
224: /* Get local mesh boundaries */
225: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
226: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
228: /* Set local loop dimensions */
229: xe = xs+xm;
230: ye = ys+ym;
231: if (xs == 0) xsm = xs-1;
232: else xsm = xs;
233: if (ys == 0) ysm = ys-1;
234: else ysm = ys;
235: if (xe == mx) xep = xe+1;
236: else xep = xe;
237: if (ye == my) yep = ye+1;
238: else yep = ye;
240: /* Compute local gradient contributions over the lower triangular elements */
241: for (j=ysm; j<ye; j++) { /* for (j=-1; j<my; j++) */
242: for (i=xsm; i<xe; i++) { /* for (i=-1; i<mx; i++) */
243: k = (j-gys)*gxm + i-gxs;
244: v = zero;
245: vr = zero;
246: vt = zero;
247: if (i >= 0 && j >= 0) v = x[k];
248: if (i < mx-1 && j > -1) vr = x[k+1];
249: if (i > -1 && j < my-1) vt = x[k+gxm];
250: dvdx = (vr-v)/hx;
251: dvdy = (vt-v)/hy;
252: if (i != -1 && j != -1) {
253: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
254: VecSetValuesLocal(G,1,&k,&val,ADD_VALUES);
255: }
256: if (i != mx-1 && j != -1) {
257: ind = k+1; val = dvdx/hx - cdiv3;
258: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
259: }
260: if (i != -1 && j != my-1) {
261: ind = k+gxm; val = dvdy/hy - cdiv3;
262: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
263: }
264: fquad += dvdx*dvdx + dvdy*dvdy;
265: flin -= cdiv3 * (v + vr + vt);
266: }
267: }
269: /* Compute local gradient contributions over the upper triangular elements */
270: for (j=ys; j<yep; j++) { /* for (j=0; j<=my; j++) */
271: for (i=xs; i<xep; i++) { /* for (i=0; i<=mx; i++) */
272: k = (j-gys)*gxm + i-gxs;
273: vb = zero;
274: vl = zero;
275: v = zero;
276: if (i < mx && j > 0) vb = x[k-gxm];
277: if (i > 0 && j < my) vl = x[k-1];
278: if (i < mx && j < my) v = x[k];
279: dvdx = (v-vl)/hx;
280: dvdy = (v-vb)/hy;
281: if (i != mx && j != 0) {
282: ind = k-gxm; val = - dvdy/hy - cdiv3;
283: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
284: }
285: if (i != 0 && j != my) {
286: ind = k-1; val = - dvdx/hx - cdiv3;
287: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
288: }
289: if (i != mx && j != my) {
290: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
291: VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
292: }
293: fquad += dvdx*dvdx + dvdy*dvdy;
294: flin -= cdiv3 * (vb + vl + v);
295: }
296: }
298: /* Restore vector */
299: VecRestoreArray(localX,&x);
301: /* Assemble gradient vector */
302: VecAssemblyBegin(G);
303: VecAssemblyEnd(G);
305: /* Scale the gradient */
306: area = p5*hx*hy;
307: floc = area * (p5 * fquad + flin);
308: VecScale(G, area);
310: /* Sum function contributions from all processes */ /* TODO: Change to PetscCallMPI() */
311: (PetscErrorCode)MPI_Allreduce((void*)&floc,(void*)f,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
313: PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16);
314: return 0;
315: }
317: PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void*ctx)
318: {
319: AppCtx *user= (AppCtx*) ctx;
320: PetscInt i,j,k;
321: PetscInt col[5],row;
322: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
323: PetscReal v[5];
324: 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;
326: /* Compute the quadratic term in the objective function */
328: /*
329: Get local grid boundaries
330: */
332: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
333: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
335: for (j=ys; j<ys+ym; j++) {
337: for (i=xs; i< xs+xm; i++) {
339: row=(j-gys)*gxm + (i-gxs);
341: k=0;
342: if (j>gys) {
343: v[k]=-2*hyhy; col[k]=row - gxm; k++;
344: }
346: if (i>gxs) {
347: v[k]= -2*hxhx; col[k]=row - 1; k++;
348: }
350: v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;
352: if (i+1 < gxs+gxm) {
353: v[k]= -2.0*hxhx; col[k]=row+1; k++;
354: }
356: if (j+1 <gys+gym) {
357: v[k]= -2*hyhy; col[k] = row+gxm; k++;
358: }
360: MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES);
362: }
363: }
364: /*
365: Assemble matrix, using the 2-step process:
366: MatAssemblyBegin(), MatAssemblyEnd().
367: By placing code between these two statements, computations can be
368: done while messages are in transition.
369: */
370: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
371: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
372: /*
373: Tell the matrix we will never add a new nonzero location to the
374: matrix. If we do it will generate an error.
375: */
376: MatScale(A,area);
377: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
378: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
379: PetscLogFlops(9*xm*ym+49*xm);
380: return 0;
381: }
383: /*TEST
385: build:
386: requires: !complex
388: test:
389: suffix: 1
390: nsize: 2
391: args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4
393: TEST*/