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: } Ctx;
13: int main(int argc, char **argv)
14: {
15: TS ts; /* time integration context */
16: Vec X; /* solution, residual vectors */
17: Mat J; /* Jacobian matrix */
18: PetscScalar *x;
19: PetscReal ftime;
20: PetscInt i, steps, nits, lits;
21: PetscBool view_final;
22: Ctx ctx;
24: PetscFunctionBeginUser;
25: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
26: ctx.n = 3;
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &ctx.n, NULL));
28: PetscCheck(ctx.n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The dimension specified with -n must be at least 2");
30: view_final = PETSC_FALSE;
31: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_final", &view_final, NULL));
33: /*
34: Create Jacobian matrix data structure and state vector
35: */
36: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
37: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx.n, ctx.n));
38: PetscCall(MatSetFromOptions(J));
39: PetscCall(MatSetUp(J));
40: PetscCall(MatCreateVecs(J, &X, NULL));
42: /* Create time integration context */
43: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
44: PetscCall(TSSetType(ts, TSPSEUDO));
45: PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &ctx));
46: PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &ctx));
47: PetscCall(TSSetMaxSteps(ts, 1000));
48: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
49: PetscCall(TSSetTimeStep(ts, 1e-3));
50: PetscCall(TSMonitorSet(ts, MonitorObjective, &ctx, NULL));
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Customize time integrator; set runtime options
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: PetscCall(TSSetFromOptions(ts));
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Evaluate initial guess; then solve nonlinear system
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: PetscCall(VecSet(X, 0.0));
61: PetscCall(VecGetArray(X, &x));
62: #if 1
63: x[0] = 5.;
64: x[1] = -5.;
65: for (i = 2; i < ctx.n; i++) x[i] = 5.;
66: #else
67: x[0] = 1.0;
68: x[1] = 15.0;
69: for (i = 2; i < ctx.n; i++) x[i] = 10.0;
70: #endif
71: PetscCall(VecRestoreArray(X, &x));
73: PetscCall(TSSolve(ts, X));
74: PetscCall(TSGetSolveTime(ts, &ftime));
75: PetscCall(TSGetStepNumber(ts, &steps));
76: PetscCall(TSGetSNESIterations(ts, &nits));
77: PetscCall(TSGetKSPIterations(ts, &lits));
78: 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));
79: if (view_final) PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Free work space. All PETSc objects should be destroyed when they
83: are no longer needed.
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: PetscCall(VecDestroy(&X));
87: PetscCall(MatDestroy(&J));
88: PetscCall(TSDestroy(&ts));
89: PetscCall(PetscFinalize());
90: return 0;
91: }
93: static PetscErrorCode MonitorObjective(TS ts, PetscInt step, PetscReal t, Vec X, void *ictx)
94: {
95: Ctx *ctx = (Ctx *)ictx;
96: const PetscScalar *x;
97: PetscScalar f;
98: PetscReal dt, gnorm;
99: PetscInt i, snesit, linit;
100: SNES snes;
101: Vec Xdot, F;
103: PetscFunctionBeginUser;
104: /* Compute objective functional */
105: PetscCall(VecGetArrayRead(X, &x));
106: f = 0;
107: for (i = 0; i < ctx->n - 1; i++) f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i + 1] - PetscSqr(x[i]));
108: PetscCall(VecRestoreArrayRead(X, &x));
110: /* Compute norm of gradient */
111: PetscCall(VecDuplicate(X, &Xdot));
112: PetscCall(VecDuplicate(X, &F));
113: PetscCall(FormIFunction(ts, t, X, Xdot, F, ictx));
114: PetscCall(VecNorm(F, NORM_2, &gnorm));
115: PetscCall(VecDestroy(&Xdot));
116: PetscCall(VecDestroy(&F));
118: PetscCall(TSGetTimeStep(ts, &dt));
119: PetscCall(TSGetSNES(ts, &snes));
120: PetscCall(SNESGetIterationNumber(snes, &snesit));
121: PetscCall(SNESGetLinearSolveIterations(snes, &linit));
122: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%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));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /* ------------------------------------------------------------------- */
127: /*
128: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
130: Input Parameters:
131: + ts - the TS context
132: . t - time
133: . X - input vector
134: . Xdot - time derivative
135: - ctx - optional user-defined context
137: Output Parameter:
138: . F - function vector
139: */
140: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ictx)
141: {
142: const PetscScalar *x;
143: PetscScalar *f;
144: Ctx *ctx = (Ctx *)ictx;
146: PetscFunctionBeginUser;
147: /*
148: Get pointers to vector data.
149: - For default PETSc vectors, VecGetArray() returns a pointer to
150: the data array. Otherwise, the routine is implementation dependent.
151: - You MUST call VecRestoreArray() when you no longer need access to
152: the array.
153: */
154: PetscCall(VecGetArrayRead(X, &x));
155: PetscCall(VecZeroEntries(F));
156: PetscCall(VecGetArray(F, &f));
158: /* Compute gradient of objective */
159: for (PetscInt i = 0; i < ctx->n - 1; i++) {
160: PetscScalar a, a0, a1;
161: a = x[i + 1] - PetscSqr(x[i]);
162: a0 = -2. * x[i];
163: a1 = 1.;
164: f[i] += -2. * (1. - x[i]) + 200. * a * a0;
165: f[i + 1] += 200. * a * a1;
166: }
167: /* Restore vectors */
168: PetscCall(VecRestoreArrayRead(X, &x));
169: PetscCall(VecRestoreArray(F, &f));
170: PetscCall(VecAXPY(F, 1.0, Xdot));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
173: /* ------------------------------------------------------------------- */
174: /*
175: FormIJacobian - Evaluates Jacobian matrix.
177: Input Parameters:
178: + ts - the TS context
179: . t - pseudo-time
180: . X - input vector
181: . Xdot - time derivative
182: . shift - multiplier for mass matrix
183: . dummy - user-defined context
185: Output Parameters:
186: . J - Jacobian matrix
187: . B - optionally different matrix used to construct the preconditioner
188: */
189: static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx)
190: {
191: const PetscScalar *x;
192: Ctx *ctx = (Ctx *)ictx;
194: PetscFunctionBeginUser;
195: PetscCall(MatZeroEntries(B));
196: /*
197: Get pointer to vector data
198: */
199: PetscCall(VecGetArrayRead(X, &x));
201: /*
202: Compute Jacobian entries and insert into matrix.
203: */
204: for (PetscInt i = 0; i < ctx->n - 1; i++) {
205: PetscInt rowcol[2];
206: PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11;
207: rowcol[0] = i;
208: rowcol[1] = i + 1;
209: a = x[i + 1] - PetscSqr(x[i]);
210: a0 = -2. * x[i];
211: a00 = -2.;
212: a01 = 0.;
213: a1 = 1.;
214: a10 = 0.;
215: a11 = 0.;
216: v[0][0] = 2. + 200. * (a * a00 + a0 * a0);
217: v[0][1] = 200. * (a * a01 + a1 * a0);
218: v[1][0] = 200. * (a * a10 + a0 * a1);
219: v[1][1] = 200. * (a * a11 + a1 * a1);
220: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES));
221: }
222: for (PetscInt i = 0; i < ctx->n; i++) PetscCall(MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES));
224: PetscCall(VecRestoreArrayRead(X, &x));
226: /*
227: Assemble matrix
228: */
229: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
230: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
231: if (J != B) {
232: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
233: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*TEST
240: test:
241: requires: !single
243: test:
244: args: -pc_type lu -ts_time_step 1e-5 -ts_max_time 1e5 -n 50 -snes_max_it 5 -snes_type newtonls -ts_max_snes_failures unlimited
245: requires: !single
246: suffix: 2
248: TEST*/