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;
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
27: ctx.n = 3;
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &ctx.n, NULL));
29: PetscCheck(ctx.n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The dimension specified with -n must be at least 2");
31: view_final = PETSC_FALSE;
32: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_final", &view_final, NULL));
34: ctx.monitor_short = PETSC_FALSE;
35: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor_short", &ctx.monitor_short, NULL));
37: /*
38: Create Jacobian matrix data structure and state vector
39: */
40: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
41: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx.n, ctx.n));
42: PetscCall(MatSetFromOptions(J));
43: PetscCall(MatSetUp(J));
44: PetscCall(MatCreateVecs(J, &X, NULL));
46: /* Create time integration context */
47: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
48: PetscCall(TSSetType(ts, TSPSEUDO));
49: PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &ctx));
50: PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &ctx));
51: PetscCall(TSSetMaxSteps(ts, 1000));
52: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
53: PetscCall(TSSetTimeStep(ts, 1e-3));
54: PetscCall(TSMonitorSet(ts, MonitorObjective, &ctx, NULL));
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Customize time integrator; set runtime options
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: PetscCall(TSSetFromOptions(ts));
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Evaluate initial guess; then solve nonlinear system
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscCall(VecSet(X, 0.0));
65: PetscCall(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: PetscCall(VecRestoreArray(X, &x));
77: PetscCall(TSSolve(ts, X));
78: PetscCall(TSGetSolveTime(ts, &ftime));
79: PetscCall(TSGetStepNumber(ts, &steps));
80: PetscCall(TSGetSNESIterations(ts, &nits));
81: PetscCall(TSGetKSPIterations(ts, &lits));
82: PetscCall(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) PetscCall(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: PetscCall(VecDestroy(&X));
91: PetscCall(MatDestroy(&J));
92: PetscCall(TSDestroy(&ts));
93: PetscCall(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;
107: PetscFunctionBeginUser;
108: /* Compute objective functional */
109: PetscCall(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: PetscCall(VecRestoreArrayRead(X, &x));
114: /* Compute norm of gradient */
115: PetscCall(VecDuplicate(X, &Xdot));
116: PetscCall(VecDuplicate(X, &F));
117: PetscCall(VecZeroEntries(Xdot));
118: PetscCall(FormIFunction(ts, t, X, Xdot, F, ictx));
119: PetscCall(VecNorm(F, NORM_2, &gnorm));
120: PetscCall(VecDestroy(&Xdot));
121: PetscCall(VecDestroy(&F));
123: PetscCall(TSGetTimeStep(ts, &dt));
124: PetscCall(TSGetSNES(ts, &snes));
125: PetscCall(SNESGetIterationNumber(snes, &snesit));
126: PetscCall(SNESGetLinearSolveIterations(snes, &linit));
127: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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;
152: PetscFunctionBeginUser;
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: PetscCall(VecGetArrayRead(X, &x));
161: PetscCall(VecZeroEntries(F));
162: PetscCall(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: PetscCall(VecRestoreArrayRead(X, &x));
175: PetscCall(VecRestoreArray(F, &f));
176: PetscCall(VecAXPY(F, 1.0, Xdot));
177: PetscFunctionReturn(PETSC_SUCCESS);
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;
202: PetscFunctionBeginUser;
203: PetscCall(MatZeroEntries(B));
204: /*
205: Get pointer to vector data
206: */
207: PetscCall(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: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES));
229: }
230: for (i = 0; i < ctx->n; i++) PetscCall(MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES));
232: PetscCall(VecRestoreArrayRead(X, &x));
234: /*
235: Assemble matrix
236: */
237: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
238: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
239: if (J != B) {
240: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
241: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
242: }
243: PetscFunctionReturn(PETSC_SUCCESS);
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*/