1 static char help[] = "Pseudotransient continuation to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";
2
3 #include <petscts.h>
4
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 *);
8
9 typedef struct {
10 PetscInt n;
11 PetscBool monitor_short;
12 } Ctx;
13
main(int argc,char ** argv)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;
24
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");
30
31 view_final = PETSC_FALSE;
32 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_final", &view_final, NULL));
33
34 ctx.monitor_short = PETSC_FALSE;
35 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor_short", &ctx.monitor_short, NULL));
36
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));
45
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));
55
56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 Customize time integrator; set runtime options
58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 PetscCall(TSSetFromOptions(ts));
60
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));
76
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));
84
85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 Free work space. All PETSc objects should be destroyed when they
87 are no longer needed.
88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89
90 PetscCall(VecDestroy(&X));
91 PetscCall(MatDestroy(&J));
92 PetscCall(TSDestroy(&ts));
93 PetscCall(PetscFinalize());
94 return 0;
95 }
96
MonitorObjective(TS ts,PetscInt step,PetscReal t,Vec X,void * ictx)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;
106
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));
113
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));
122
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 }
130
131 /* ------------------------------------------------------------------- */
132 /*
133 FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
134
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
141
142 Output Parameter:
143 . F - function vector
144 */
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ictx)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;
151
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));
163
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.
182
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
190
191 Output Parameters:
192 . J - Jacobian matrix
193 . B - optionally different matrix used to construct the preconditioner
194 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat J,Mat B,void * ictx)195 static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx)
196 {
197 const PetscScalar *x;
198 PetscInt i;
199 Ctx *ctx = (Ctx *)ictx;
200
201 PetscFunctionBeginUser;
202 PetscCall(MatZeroEntries(B));
203 /*
204 Get pointer to vector data
205 */
206 PetscCall(VecGetArrayRead(X, &x));
207
208 /*
209 Compute Jacobian entries and insert into matrix.
210 */
211 for (i = 0; i < ctx->n - 1; i++) {
212 PetscInt rowcol[2];
213 PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11;
214 rowcol[0] = i;
215 rowcol[1] = i + 1;
216 a = x[i + 1] - PetscSqr(x[i]);
217 a0 = -2. * x[i];
218 a00 = -2.;
219 a01 = 0.;
220 a1 = 1.;
221 a10 = 0.;
222 a11 = 0.;
223 v[0][0] = 2. + 200. * (a * a00 + a0 * a0);
224 v[0][1] = 200. * (a * a01 + a1 * a0);
225 v[1][0] = 200. * (a * a10 + a0 * a1);
226 v[1][1] = 200. * (a * a11 + a1 * a1);
227 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES));
228 }
229 for (i = 0; i < ctx->n; i++) PetscCall(MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES));
230
231 PetscCall(VecRestoreArrayRead(X, &x));
232
233 /*
234 Assemble matrix
235 */
236 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
237 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
238 if (J != B) {
239 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
240 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
241 }
242 PetscFunctionReturn(PETSC_SUCCESS);
243 }
244
245 /*TEST
246
247 test:
248 requires: !single
249
250 test:
251 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
252 requires: !single
253 suffix: 2
254
255 TEST*/
256