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: PetscCall(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: PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg));
84: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
85: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
87: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Elastic-Plastic Torsion Problem -----\n"));
88: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
90: /* Set up distributed array */
91: PetscCall(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: PetscCall(DMSetFromOptions(user.dm));
93: PetscCall(DMSetUp(user.dm));
95: /* Create vectors */
96: PetscCall(DMCreateGlobalVector(user.dm, &x));
98: PetscCall(DMCreateLocalVector(user.dm, &user.localX));
100: /* Create Hessian */
101: PetscCall(DMCreateMatrix(user.dm, &H));
102: PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE));
104: /* The TAO code begins here */
106: /* Create TAO solver and set desired solution method */
107: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
108: PetscCall(TaoSetType(tao, TAOCG));
110: /* Set initial solution guess */
111: PetscCall(FormInitialGuess(&user, x));
112: PetscCall(TaoSetSolution(tao, x));
114: /* Set routine for function and gradient evaluation */
115: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
117: PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user));
119: /* Check for any TAO command line options */
120: PetscCall(TaoSetFromOptions(tao));
122: PetscCall(TaoGetKSP(tao, &ksp));
123: if (ksp) {
124: PetscCall(KSPGetPC(ksp, &pc));
125: PetscCall(PCSetType(pc, PCNONE));
126: }
128: /* SOLVE THE APPLICATION */
129: PetscCall(TaoSolve(tao));
131: /* Free TAO data structures */
132: PetscCall(TaoDestroy(&tao));
134: /* Free PETSc data structures */
135: PetscCall(VecDestroy(&x));
136: PetscCall(MatDestroy(&H));
138: PetscCall(VecDestroy(&user.localX));
139: PetscCall(DMDestroy(&user.dm));
141: PetscCall(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: PetscFunctionBegin;
163: /* Get local mesh boundaries */
164: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
165: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
167: /* Compute initial guess over locally owned part of mesh */
168: xe = xs + xm;
169: ye = ys + ym;
170: for (j = ys; j < ye; j++) { /* for (j=0; j<my; j++) */
171: temp = PetscMin(j + 1, my - j) * hy;
172: for (i = xs; i < xe; i++) { /* for (i=0; i<mx; i++) */
173: k = (j - gys) * gxm + i - gxs;
174: val = PetscMin((PetscMin(i + 1, mx - i)) * hx, temp);
175: PetscCall(VecSetValuesLocal(X, 1, &k, &val, ADD_VALUES));
176: }
177: }
178: PetscCall(VecAssemblyBegin(X));
179: PetscCall(VecAssemblyEnd(X));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /* ------------------------------------------------------------------ */
184: /*
185: FormFunctionGradient - Evaluates the function and corresponding gradient.
187: Input Parameters:
188: tao - the Tao context
189: X - the input vector
190: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
192: Output Parameters:
193: f - the newly evaluated function
194: G - the newly evaluated gradient
195: */
196: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
197: {
198: AppCtx *user = (AppCtx *)ptr;
199: PetscInt i, j, k, ind;
200: PetscInt xe, ye, xsm, ysm, xep, yep;
201: PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
202: PetscInt mx = user->mx, my = user->my;
203: PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param / three;
204: PetscReal p5 = 0.5, area, val, flin, fquad;
205: PetscReal v, vb, vl, vr, vt, dvdx, dvdy;
206: PetscReal hx = 1.0 / (user->mx + 1);
207: PetscReal hy = 1.0 / (user->my + 1);
208: Vec localX = user->localX;
210: PetscFunctionBegin;
211: /* Initialize */
212: flin = fquad = zero;
214: PetscCall(VecSet(G, zero));
215: /*
216: Scatter ghost points to local vector,using the 2-step process
217: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
218: By placing code between these two statements, computations can be
219: done while messages are in transition.
220: */
221: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
222: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
224: /* Get pointer to vector data */
225: PetscCall(VecGetArray(localX, &x));
227: /* Get local mesh boundaries */
228: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
229: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
231: /* Set local loop dimensions */
232: xe = xs + xm;
233: ye = ys + ym;
234: if (xs == 0) xsm = xs - 1;
235: else xsm = xs;
236: if (ys == 0) ysm = ys - 1;
237: else ysm = ys;
238: if (xe == mx) xep = xe + 1;
239: else xep = xe;
240: if (ye == my) yep = ye + 1;
241: else yep = ye;
243: /* Compute local gradient contributions over the lower triangular elements */
244: for (j = ysm; j < ye; j++) { /* for (j=-1; j<my; j++) */
245: for (i = xsm; i < xe; i++) { /* for (i=-1; i<mx; i++) */
246: k = (j - gys) * gxm + i - gxs;
247: v = zero;
248: vr = zero;
249: vt = zero;
250: if (i >= 0 && j >= 0) v = x[k];
251: if (i < mx - 1 && j > -1) vr = x[k + 1];
252: if (i > -1 && j < my - 1) vt = x[k + gxm];
253: dvdx = (vr - v) / hx;
254: dvdy = (vt - v) / hy;
255: if (i != -1 && j != -1) {
256: ind = k;
257: val = -dvdx / hx - dvdy / hy - cdiv3;
258: PetscCall(VecSetValuesLocal(G, 1, &k, &val, ADD_VALUES));
259: }
260: if (i != mx - 1 && j != -1) {
261: ind = k + 1;
262: val = dvdx / hx - cdiv3;
263: PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES));
264: }
265: if (i != -1 && j != my - 1) {
266: ind = k + gxm;
267: val = dvdy / hy - cdiv3;
268: PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES));
269: }
270: fquad += dvdx * dvdx + dvdy * dvdy;
271: flin -= cdiv3 * (v + vr + vt);
272: }
273: }
275: /* Compute local gradient contributions over the upper triangular elements */
276: for (j = ys; j < yep; j++) { /* for (j=0; j<=my; j++) */
277: for (i = xs; i < xep; i++) { /* for (i=0; i<=mx; i++) */
278: k = (j - gys) * gxm + i - gxs;
279: vb = zero;
280: vl = zero;
281: v = zero;
282: if (i < mx && j > 0) vb = x[k - gxm];
283: if (i > 0 && j < my) vl = x[k - 1];
284: if (i < mx && j < my) v = x[k];
285: dvdx = (v - vl) / hx;
286: dvdy = (v - vb) / hy;
287: if (i != mx && j != 0) {
288: ind = k - gxm;
289: val = -dvdy / hy - cdiv3;
290: PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES));
291: }
292: if (i != 0 && j != my) {
293: ind = k - 1;
294: val = -dvdx / hx - cdiv3;
295: PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES));
296: }
297: if (i != mx && j != my) {
298: ind = k;
299: val = dvdx / hx + dvdy / hy - cdiv3;
300: PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES));
301: }
302: fquad += dvdx * dvdx + dvdy * dvdy;
303: flin -= cdiv3 * (vb + vl + v);
304: }
305: }
307: /* Restore vector */
308: PetscCall(VecRestoreArray(localX, &x));
310: /* Assemble gradient vector */
311: PetscCall(VecAssemblyBegin(G));
312: PetscCall(VecAssemblyEnd(G));
314: /* Scale the gradient */
315: area = p5 * hx * hy;
316: floc = area * (p5 * fquad + flin);
317: PetscCall(VecScale(G, area));
319: /* Sum function contributions from all processes */ /* TODO: Change to PetscCallMPI() */
320: PetscCall((PetscErrorCode)MPI_Allreduce((void *)&floc, (void *)f, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
322: PetscCall(PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void *ctx)
327: {
328: AppCtx *user = (AppCtx *)ctx;
329: PetscInt i, j, k;
330: PetscInt col[5], row;
331: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
332: PetscReal v[5];
333: 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;
335: /* Compute the quadratic term in the objective function */
337: /*
338: Get local grid boundaries
339: */
341: PetscFunctionBegin;
342: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
343: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
345: for (j = ys; j < ys + ym; j++) {
346: for (i = xs; i < xs + xm; i++) {
347: row = (j - gys) * gxm + (i - gxs);
349: k = 0;
350: if (j > gys) {
351: v[k] = -2 * hyhy;
352: col[k] = row - gxm;
353: k++;
354: }
356: if (i > gxs) {
357: v[k] = -2 * hxhx;
358: col[k] = row - 1;
359: k++;
360: }
362: v[k] = 4.0 * (hxhx + hyhy);
363: col[k] = row;
364: k++;
366: if (i + 1 < gxs + gxm) {
367: v[k] = -2.0 * hxhx;
368: col[k] = row + 1;
369: k++;
370: }
372: if (j + 1 < gys + gym) {
373: v[k] = -2 * hyhy;
374: col[k] = row + gxm;
375: k++;
376: }
378: PetscCall(MatSetValuesLocal(A, 1, &row, k, col, v, INSERT_VALUES));
379: }
380: }
381: /*
382: Assemble matrix, using the 2-step process:
383: MatAssemblyBegin(), MatAssemblyEnd().
384: By placing code between these two statements, computations can be
385: done while messages are in transition.
386: */
387: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
388: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
389: /*
390: Tell the matrix we will never add a new nonzero location to the
391: matrix. If we do it will generate an error.
392: */
393: PetscCall(MatScale(A, area));
394: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
395: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
396: PetscCall(PetscLogFlops(9 * xm * ym + 49 * xm));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*TEST
402: build:
403: requires: !complex
405: test:
406: suffix: 1
407: nsize: 2
408: args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4
410: TEST*/