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, NULL, 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(FormIFunction(ts, t, X, Xdot, F, ictx));
118: PetscCall(VecNorm(F, NORM_2, &gnorm));
119: PetscCall(VecDestroy(&Xdot));
120: PetscCall(VecDestroy(&F));
122: PetscCall(TSGetTimeStep(ts, &dt));
123: PetscCall(TSGetSNES(ts, &snes));
124: PetscCall(SNESGetIterationNumber(snes, &snesit));
125: PetscCall(SNESGetLinearSolveIterations(snes, &linit));
126: 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));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /* ------------------------------------------------------------------- */
131: /*
132: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
134: Input Parameters:
135: + ts - the TS context
136: . t - time
137: . X - input vector
138: . Xdot - time derivative
139: - ctx - optional user-defined context
141: Output Parameter:
142: . F - function vector
143: */
144: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ictx)
145: {
146: const PetscScalar *x;
147: PetscScalar *f;
148: Ctx *ctx = (Ctx *)ictx;
150: PetscFunctionBeginUser;
151: /*
152: Get pointers to vector data.
153: - For default PETSc vectors, VecGetArray() returns a pointer to
154: the data array. Otherwise, the routine is implementation dependent.
155: - You MUST call VecRestoreArray() when you no longer need access to
156: the array.
157: */
158: PetscCall(VecGetArrayRead(X, &x));
159: PetscCall(VecZeroEntries(F));
160: PetscCall(VecGetArray(F, &f));
162: /* Compute gradient of objective */
163: for (PetscInt i = 0; i < ctx->n - 1; i++) {
164: PetscScalar a, a0, a1;
165: a = x[i + 1] - PetscSqr(x[i]);
166: a0 = -2. * x[i];
167: a1 = 1.;
168: f[i] += -2. * (1. - x[i]) + 200. * a * a0;
169: f[i + 1] += 200. * a * a1;
170: }
171: /* Restore vectors */
172: PetscCall(VecRestoreArrayRead(X, &x));
173: PetscCall(VecRestoreArray(F, &f));
174: PetscCall(VecAXPY(F, 1.0, Xdot));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
177: /* ------------------------------------------------------------------- */
178: /*
179: FormIJacobian - Evaluates Jacobian matrix.
181: Input Parameters:
182: + ts - the TS context
183: . t - pseudo-time
184: . X - input vector
185: . Xdot - time derivative
186: . shift - multiplier for mass matrix
187: . dummy - user-defined context
189: Output Parameters:
190: . J - Jacobian matrix
191: . B - optionally different matrix used to construct the preconditioner
192: */
193: static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx)
194: {
195: const PetscScalar *x;
196: Ctx *ctx = (Ctx *)ictx;
198: PetscFunctionBeginUser;
199: PetscCall(MatZeroEntries(B));
200: /*
201: Get pointer to vector data
202: */
203: PetscCall(VecGetArrayRead(X, &x));
205: /*
206: Compute Jacobian entries and insert into matrix.
207: */
208: for (PetscInt i = 0; i < ctx->n - 1; i++) {
209: PetscInt rowcol[2];
210: PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11;
211: rowcol[0] = i;
212: rowcol[1] = i + 1;
213: a = x[i + 1] - PetscSqr(x[i]);
214: a0 = -2. * x[i];
215: a00 = -2.;
216: a01 = 0.;
217: a1 = 1.;
218: a10 = 0.;
219: a11 = 0.;
220: v[0][0] = 2. + 200. * (a * a00 + a0 * a0);
221: v[0][1] = 200. * (a * a01 + a1 * a0);
222: v[1][0] = 200. * (a * a10 + a0 * a1);
223: v[1][1] = 200. * (a * a11 + a1 * a1);
224: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES));
225: }
226: for (PetscInt i = 0; i < ctx->n; i++) PetscCall(MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES));
228: PetscCall(VecRestoreArrayRead(X, &x));
230: /*
231: Assemble matrix
232: */
233: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
234: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
235: if (J != B) {
236: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
237: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
238: }
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*TEST
244: test:
245: requires: !single
247: test:
248: args: -pc_type lu -ts_time_step 1e-5 -ts_max_time 1e5 -n 50 -monitor_short -snes_max_it 5 -snes_type newtonls -ts_max_snes_failures unlimited
249: requires: !single
250: suffix: 2
252: TEST*/