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[] = "Demonstrates use of the TAO package to solve \n\
37: unconstrained minimization problems in parallel. This example is based on \n\
38: the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.\n\
39: The command line options are:\n\
40: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
41: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
42: -par <param>, where <param> = angle of twist per unit length\n\n";
44: /*
45: User-defined application context - contains data needed by the
46: application-provided call-back routines, FormFunction() and
47: FormGradient().
48: */
49: typedef struct {
50: /* parameters */
51: PetscInt mx, my; /* global discretization in x- and y-directions */
52: PetscReal param; /* nonlinearity parameter */
54: /* work space */
55: Vec localX; /* local vectors */
56: DM dm; /* distributed array data structure */
57: } AppCtx;
59: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
60: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
61: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
63: int main(int argc, char **argv)
64: {
65: Vec x;
66: Mat H;
67: PetscInt Nx, Ny;
68: Tao tao;
69: PetscBool flg;
70: KSP ksp;
71: PC pc;
72: AppCtx user;
74: PetscInitialize(&argc, &argv, (char *)0, help);
76: /* Specify default dimension of the problem */
77: user.param = 5.0;
78: user.mx = 10;
79: user.my = 10;
80: Nx = Ny = PETSC_DECIDE;
82: /* Check for any command line arguments that override defaults */
83: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg);
84: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);
85: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
87: PetscPrintf(PETSC_COMM_WORLD, "\n---- Elastic-Plastic Torsion Problem -----\n");
88: PetscPrintf(PETSC_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my);
90: /* Set up distributed array */
91: 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);
92: DMSetFromOptions(user.dm);
93: DMSetUp(user.dm);
95: /* Create vectors */
96: DMCreateGlobalVector(user.dm, &x);
98: DMCreateLocalVector(user.dm, &user.localX);
100: /* Create Hessian */
101: DMCreateMatrix(user.dm, &H);
102: MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);
104: /* The TAO code begins here */
106: /* Create TAO solver and set desired solution method */
107: TaoCreate(PETSC_COMM_WORLD, &tao);
108: TaoSetType(tao, TAOCG);
110: /* Set initial solution guess */
111: FormInitialGuess(&user, x);
112: TaoSetSolution(tao, x);
114: /* Set routine for function and gradient evaluation */
115: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
117: TaoSetHessian(tao, H, H, FormHessian, (void *)&user);
119: /* Check for any TAO command line options */
120: TaoSetFromOptions(tao);
122: TaoGetKSP(tao, &ksp);
123: if (ksp) {
124: KSPGetPC(ksp, &pc);
125: PCSetType(pc, PCNONE);
126: }
128: /* SOLVE THE APPLICATION */
129: TaoSolve(tao);
131: /* Free TAO data structures */
132: TaoDestroy(&tao);
134: /* Free PETSc data structures */
135: VecDestroy(&x);
136: MatDestroy(&H);
138: VecDestroy(&user.localX);
139: DMDestroy(&user.dm);
141: PetscFinalize();
142: return 0;
143: }
145: /* ------------------------------------------------------------------- */
146: /*
147: FormInitialGuess - Computes an initial approximation to the solution.
149: Input Parameters:
150: . user - user-defined application context
151: . X - vector
153: Output Parameters:
154: X - vector
155: */
156: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
157: {
158: PetscInt i, j, k, mx = user->mx, my = user->my;
159: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
160: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), temp, val;
162: /* Get local mesh boundaries */
163: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
164: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
166: /* Compute initial guess over locally owned part of mesh */
167: xe = xs + xm;
168: ye = ys + ym;
169: for (j = ys; j < ye; j++) { /* for (j=0; j<my; j++) */
170: temp = PetscMin(j + 1, my - j) * hy;
171: for (i = xs; i < xe; i++) { /* for (i=0; i<mx; i++) */
172: k = (j - gys) * gxm + i - gxs;
173: val = PetscMin((PetscMin(i + 1, mx - i)) * hx, temp);
174: VecSetValuesLocal(X, 1, &k, &val, ADD_VALUES);
175: }
176: }
177: VecAssemblyBegin(X);
178: VecAssemblyEnd(X);
179: return 0;
180: }
182: /* ------------------------------------------------------------------ */
183: /*
184: FormFunctionGradient - Evaluates the function and corresponding gradient.
186: Input Parameters:
187: tao - the Tao context
188: X - the input vector
189: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
191: Output Parameters:
192: f - the newly evaluated function
193: G - the newly evaluated gradient
194: */
195: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
196: {
197: AppCtx *user = (AppCtx *)ptr;
198: PetscInt i, j, k, ind;
199: PetscInt xe, ye, xsm, ysm, xep, yep;
200: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
201: PetscInt mx = user->mx, my = user->my;
202: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param / three;
203: PetscReal p5 = 0.5, area, val, flin, fquad;
204: PetscReal v, vb, vl, vr, vt, dvdx, dvdy;
205: PetscReal hx = 1.0 / (user->mx + 1);
206: PetscReal hy = 1.0 / (user->my + 1);
207: Vec localX = user->localX;
209: /* Initialize */
210: flin = fquad = zero;
212: VecSet(G, zero);
213: /*
214: Scatter ghost points to local vector,using the 2-step process
215: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
216: By placing code between these two statements, computations can be
217: done while messages are in transition.
218: */
219: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
220: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
222: /* Get pointer to vector data */
223: VecGetArray(localX, &x);
225: /* Get local mesh boundaries */
226: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
227: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
229: /* Set local loop dimensions */
230: xe = xs + xm;
231: ye = ys + ym;
232: if (xs == 0) xsm = xs - 1;
233: else xsm = xs;
234: if (ys == 0) ysm = ys - 1;
235: else ysm = ys;
236: if (xe == mx) xep = xe + 1;
237: else xep = xe;
238: if (ye == my) yep = ye + 1;
239: else yep = ye;
241: /* Compute local gradient contributions over the lower triangular elements */
242: for (j = ysm; j < ye; j++) { /* for (j=-1; j<my; j++) */
243: for (i = xsm; i < xe; i++) { /* for (i=-1; i<mx; i++) */
244: k = (j - gys) * gxm + i - gxs;
245: v = zero;
246: vr = zero;
247: vt = zero;
248: if (i >= 0 && j >= 0) v = x[k];
249: if (i < mx - 1 && j > -1) vr = x[k + 1];
250: if (i > -1 && j < my - 1) vt = x[k + gxm];
251: dvdx = (vr - v) / hx;
252: dvdy = (vt - v) / hy;
253: if (i != -1 && j != -1) {
254: ind = k;
255: val = -dvdx / hx - dvdy / hy - cdiv3;
256: VecSetValuesLocal(G, 1, &k, &val, ADD_VALUES);
257: }
258: if (i != mx - 1 && j != -1) {
259: ind = k + 1;
260: val = dvdx / hx - cdiv3;
261: VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES);
262: }
263: if (i != -1 && j != my - 1) {
264: ind = k + gxm;
265: val = dvdy / hy - cdiv3;
266: VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES);
267: }
268: fquad += dvdx * dvdx + dvdy * dvdy;
269: flin -= cdiv3 * (v + vr + vt);
270: }
271: }
273: /* Compute local gradient contributions over the upper triangular elements */
274: for (j = ys; j < yep; j++) { /* for (j=0; j<=my; j++) */
275: for (i = xs; i < xep; i++) { /* for (i=0; i<=mx; i++) */
276: k = (j - gys) * gxm + i - gxs;
277: vb = zero;
278: vl = zero;
279: v = zero;
280: if (i < mx && j > 0) vb = x[k - gxm];
281: if (i > 0 && j < my) vl = x[k - 1];
282: if (i < mx && j < my) v = x[k];
283: dvdx = (v - vl) / hx;
284: dvdy = (v - vb) / hy;
285: if (i != mx && j != 0) {
286: ind = k - gxm;
287: val = -dvdy / hy - cdiv3;
288: VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES);
289: }
290: if (i != 0 && j != my) {
291: ind = k - 1;
292: val = -dvdx / hx - cdiv3;
293: VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES);
294: }
295: if (i != mx && j != my) {
296: ind = k;
297: val = dvdx / hx + dvdy / hy - cdiv3;
298: VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES);
299: }
300: fquad += dvdx * dvdx + dvdy * dvdy;
301: flin -= cdiv3 * (vb + vl + v);
302: }
303: }
305: /* Restore vector */
306: VecRestoreArray(localX, &x);
308: /* Assemble gradient vector */
309: VecAssemblyBegin(G);
310: VecAssemblyEnd(G);
312: /* Scale the gradient */
313: area = p5 * hx * hy;
314: floc = area * (p5 * fquad + flin);
315: VecScale(G, area);
317: /* Sum function contributions from all processes */ /* TODO: Change to PetscCallMPI() */
318: (PetscErrorCode)MPI_Allreduce((void *)&floc, (void *)f, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD);
320: PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16);
321: return 0;
322: }
324: PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void *ctx)
325: {
326: AppCtx *user = (AppCtx *)ctx;
327: PetscInt i, j, k;
328: PetscInt col[5], row;
329: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
330: PetscReal v[5];
331: 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;
333: /* Compute the quadratic term in the objective function */
335: /*
336: Get local grid boundaries
337: */
339: DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
340: DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
342: for (j = ys; j < ys + ym; j++) {
343: for (i = xs; i < xs + xm; i++) {
344: row = (j - gys) * gxm + (i - gxs);
346: k = 0;
347: if (j > gys) {
348: v[k] = -2 * hyhy;
349: col[k] = row - gxm;
350: k++;
351: }
353: if (i > gxs) {
354: v[k] = -2 * hxhx;
355: col[k] = row - 1;
356: k++;
357: }
359: v[k] = 4.0 * (hxhx + hyhy);
360: col[k] = row;
361: k++;
363: if (i + 1 < gxs + gxm) {
364: v[k] = -2.0 * hxhx;
365: col[k] = row + 1;
366: k++;
367: }
369: if (j + 1 < gys + gym) {
370: v[k] = -2 * hyhy;
371: col[k] = row + gxm;
372: k++;
373: }
375: MatSetValuesLocal(A, 1, &row, k, col, v, INSERT_VALUES);
376: }
377: }
378: /*
379: Assemble matrix, using the 2-step process:
380: MatAssemblyBegin(), MatAssemblyEnd().
381: By placing code between these two statements, computations can be
382: done while messages are in transition.
383: */
384: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
385: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
386: /*
387: Tell the matrix we will never add a new nonzero location to the
388: matrix. If we do it will generate an error.
389: */
390: MatScale(A, area);
391: MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
392: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
393: PetscLogFlops(9 * xm * ym + 49 * xm);
394: return 0;
395: }
397: /*TEST
399: build:
400: requires: !complex
402: test:
403: suffix: 1
404: nsize: 2
405: args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4
407: TEST*/