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