Actual source code: ex24.c
1: static char help[] = "Pseudotransient continuation to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";
3: #include <petscts.h>
5: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
6: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
7: static PetscErrorCode MonitorObjective(TS, PetscInt, PetscReal, Vec, void *);
9: typedef struct {
10: PetscInt n;
11: PetscBool monitor_short;
12: } Ctx;
14: int main(int argc, char **argv)
15: {
16: TS ts; /* time integration context */
17: Vec X; /* solution, residual vectors */
18: Mat J; /* Jacobian matrix */
19: PetscScalar *x;
20: PetscReal ftime;
21: PetscInt i, steps, nits, lits;
22: PetscBool view_final;
23: Ctx ctx;
26: PetscInitialize(&argc, &argv, (char *)0, help);
27: ctx.n = 3;
28: PetscOptionsGetInt(NULL, NULL, "-n", &ctx.n, NULL);
31: view_final = PETSC_FALSE;
32: PetscOptionsGetBool(NULL, NULL, "-view_final", &view_final, NULL);
34: ctx.monitor_short = PETSC_FALSE;
35: PetscOptionsGetBool(NULL, NULL, "-monitor_short", &ctx.monitor_short, NULL);
37: /*
38: Create Jacobian matrix data structure and state vector
39: */
40: MatCreate(PETSC_COMM_WORLD, &J);
41: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx.n, ctx.n);
42: MatSetFromOptions(J);
43: MatSetUp(J);
44: MatCreateVecs(J, &X, NULL);
46: /* Create time integration context */
47: TSCreate(PETSC_COMM_WORLD, &ts);
48: TSSetType(ts, TSPSEUDO);
49: TSSetIFunction(ts, NULL, FormIFunction, &ctx);
50: TSSetIJacobian(ts, J, J, FormIJacobian, &ctx);
51: TSSetMaxSteps(ts, 1000);
52: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
53: TSSetTimeStep(ts, 1e-3);
54: TSMonitorSet(ts, MonitorObjective, &ctx, NULL);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Customize time integrator; set runtime options
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: TSSetFromOptions(ts);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Evaluate initial guess; then solve nonlinear system
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: VecSet(X, 0.0);
65: VecGetArray(X, &x);
66: #if 1
67: x[0] = 5.;
68: x[1] = -5.;
69: for (i = 2; i < ctx.n; i++) x[i] = 5.;
70: #else
71: x[0] = 1.0;
72: x[1] = 15.0;
73: for (i = 2; i < ctx.n; i++) x[i] = 10.0;
74: #endif
75: VecRestoreArray(X, &x);
77: TSSolve(ts, X);
78: TSGetSolveTime(ts, &ftime);
79: TSGetStepNumber(ts, &steps);
80: TSGetSNESIterations(ts, &nits);
81: TSGetKSPIterations(ts, &lits);
82: PetscPrintf(PETSC_COMM_WORLD, "Time integrator took (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") iterations to reach final time %g\n", steps, nits, lits, (double)ftime);
83: if (view_final) VecView(X, PETSC_VIEWER_STDOUT_WORLD);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Free work space. All PETSc objects should be destroyed when they
87: are no longer needed.
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: VecDestroy(&X);
91: MatDestroy(&J);
92: TSDestroy(&ts);
93: PetscFinalize();
94: return 0;
95: }
97: static PetscErrorCode MonitorObjective(TS ts, PetscInt step, PetscReal t, Vec X, void *ictx)
98: {
99: Ctx *ctx = (Ctx *)ictx;
100: const PetscScalar *x;
101: PetscScalar f;
102: PetscReal dt, gnorm;
103: PetscInt i, snesit, linit;
104: SNES snes;
105: Vec Xdot, F;
108: /* Compute objective functional */
109: VecGetArrayRead(X, &x);
110: f = 0;
111: for (i = 0; i < ctx->n - 1; i++) f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i + 1] - PetscSqr(x[i]));
112: VecRestoreArrayRead(X, &x);
114: /* Compute norm of gradient */
115: VecDuplicate(X, &Xdot);
116: VecDuplicate(X, &F);
117: VecZeroEntries(Xdot);
118: FormIFunction(ts, t, X, Xdot, F, ictx);
119: VecNorm(F, NORM_2, &gnorm);
120: VecDestroy(&Xdot);
121: VecDestroy(&F);
123: TSGetTimeStep(ts, &dt);
124: TSGetSNES(ts, &snes);
125: SNESGetIterationNumber(snes, &snesit);
126: SNESGetLinearSolveIterations(snes, &linit);
127: PetscPrintf(PETSC_COMM_WORLD, (ctx->monitor_short ? "%3" PetscInt_FMT " t=%10.1e dt=%10.1e f=%10.1e df=%10.1e it=(%2" PetscInt_FMT ",%3" PetscInt_FMT ")\n" : "%3" PetscInt_FMT " t=%10.4e dt=%10.4e f=%10.4e df=%10.4e it=(%2" PetscInt_FMT ",%3" PetscInt_FMT ")\n"), step, (double)t, (double)dt, (double)PetscRealPart(f), (double)gnorm, snesit, linit);
128: return 0;
129: }
131: /* ------------------------------------------------------------------- */
132: /*
133: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
135: Input Parameters:
136: + ts - the TS context
137: . t - time
138: . X - input vector
139: . Xdot - time derivative
140: - ctx - optional user-defined context
142: Output Parameter:
143: . F - function vector
144: */
145: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ictx)
146: {
147: const PetscScalar *x;
148: PetscScalar *f;
149: PetscInt i;
150: Ctx *ctx = (Ctx *)ictx;
153: /*
154: Get pointers to vector data.
155: - For default PETSc vectors, VecGetArray() returns a pointer to
156: the data array. Otherwise, the routine is implementation dependent.
157: - You MUST call VecRestoreArray() when you no longer need access to
158: the array.
159: */
160: VecGetArrayRead(X, &x);
161: VecZeroEntries(F);
162: VecGetArray(F, &f);
164: /* Compute gradient of objective */
165: for (i = 0; i < ctx->n - 1; i++) {
166: PetscScalar a, a0, a1;
167: a = x[i + 1] - PetscSqr(x[i]);
168: a0 = -2. * x[i];
169: a1 = 1.;
170: f[i] += -2. * (1. - x[i]) + 200. * a * a0;
171: f[i + 1] += 200. * a * a1;
172: }
173: /* Restore vectors */
174: VecRestoreArrayRead(X, &x);
175: VecRestoreArray(F, &f);
176: VecAXPY(F, 1.0, Xdot);
177: return 0;
178: }
179: /* ------------------------------------------------------------------- */
180: /*
181: FormIJacobian - Evaluates Jacobian matrix.
183: Input Parameters:
184: + ts - the TS context
185: . t - pseudo-time
186: . X - input vector
187: . Xdot - time derivative
188: . shift - multiplier for mass matrix
189: . dummy - user-defined context
191: Output Parameters:
192: . J - Jacobian matrix
193: . B - optionally different preconditioning matrix
194: . flag - flag indicating matrix structure
195: */
196: static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx)
197: {
198: const PetscScalar *x;
199: PetscInt i;
200: Ctx *ctx = (Ctx *)ictx;
203: MatZeroEntries(B);
204: /*
205: Get pointer to vector data
206: */
207: VecGetArrayRead(X, &x);
209: /*
210: Compute Jacobian entries and insert into matrix.
211: */
212: for (i = 0; i < ctx->n - 1; i++) {
213: PetscInt rowcol[2];
214: PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11;
215: rowcol[0] = i;
216: rowcol[1] = i + 1;
217: a = x[i + 1] - PetscSqr(x[i]);
218: a0 = -2. * x[i];
219: a00 = -2.;
220: a01 = 0.;
221: a1 = 1.;
222: a10 = 0.;
223: a11 = 0.;
224: v[0][0] = 2. + 200. * (a * a00 + a0 * a0);
225: v[0][1] = 200. * (a * a01 + a1 * a0);
226: v[1][0] = 200. * (a * a10 + a0 * a1);
227: v[1][1] = 200. * (a * a11 + a1 * a1);
228: MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES);
229: }
230: for (i = 0; i < ctx->n; i++) MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES);
232: VecRestoreArrayRead(X, &x);
234: /*
235: Assemble matrix
236: */
237: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
238: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
239: if (J != B) {
240: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
241: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
242: }
243: return 0;
244: }
246: /*TEST
248: test:
249: requires: !single
251: test:
252: args: -pc_type lu -ts_dt 1e-5 -ts_max_time 1e5 -n 50 -monitor_short -snes_max_it 5 -snes_type newtonls -ts_max_snes_failures -1
253: requires: !single
254: suffix: 2
256: TEST*/