Actual source code: ex20fwd.c
1: #define c11 1.0
2: #define c12 0
3: #define c21 2.0
4: #define c22 1.0
5: static char help[] = "Solves the van der Pol equation.\n\
6: Input parameters include:\n";
8: /* ------------------------------------------------------------------------
10: This code demonstrates how to compute trajectory sensitivities w.r.t. the stiffness parameter mu.
11: 1) Use two vectors s and sp for sensitivities w.r.t. initial values and parameters respectively. This case requires the original Jacobian matrix and a JacobianP matrix for the only parameter mu.
12: 2) Consider the initial values to be parameters as well. Then there are three parameters in total. The JacobianP matrix will be combined matrix of the Jacobian matrix and JacobianP matrix in the previous case. This choice can be selected by using command line option '-combined'
14: ------------------------------------------------------------------------- */
15: #include <petscts.h>
16: #include <petsctao.h>
18: typedef struct _n_User *User;
19: struct _n_User {
20: PetscReal mu;
21: PetscReal next_output;
22: PetscBool combined;
23: /* Sensitivity analysis support */
24: PetscInt steps;
25: PetscReal ftime;
26: Mat Jac; /* Jacobian matrix */
27: Mat Jacp; /* JacobianP matrix */
28: Vec x;
29: Mat sp; /* forward sensitivity variables */
30: };
32: /*
33: User-defined routines
34: */
35: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
36: {
37: User user = (User)ctx;
38: const PetscScalar *x, *xdot;
39: PetscScalar *f;
41: PetscFunctionBeginUser;
42: PetscCall(VecGetArrayRead(X, &x));
43: PetscCall(VecGetArrayRead(Xdot, &xdot));
44: PetscCall(VecGetArray(F, &f));
45: f[0] = xdot[0] - x[1];
46: f[1] = c21 * (xdot[0] - x[1]) + xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]);
47: PetscCall(VecRestoreArrayRead(X, &x));
48: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
49: PetscCall(VecRestoreArray(F, &f));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
54: {
55: User user = (User)ctx;
56: PetscInt rowcol[] = {0, 1};
57: PetscScalar J[2][2];
58: const PetscScalar *x;
60: PetscFunctionBeginUser;
61: PetscCall(VecGetArrayRead(X, &x));
62: J[0][0] = a;
63: J[0][1] = -1.0;
64: J[1][0] = c21 * a + user->mu * (1.0 + 2.0 * x[0] * x[1]);
65: J[1][1] = -c21 + a - user->mu * (1.0 - x[0] * x[0]);
66: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
67: PetscCall(VecRestoreArrayRead(X, &x));
69: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
71: if (A != B) {
72: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
73: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
74: }
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx)
79: {
80: User user = (User)ctx;
81: PetscInt row[] = {0, 1}, col[] = {0};
82: PetscScalar J[2][1];
83: const PetscScalar *x;
85: PetscFunctionBeginUser;
86: if (user->combined) col[0] = 2;
87: PetscCall(VecGetArrayRead(X, &x));
88: J[0][0] = 0;
89: J[1][0] = (1. - x[0] * x[0]) * x[1] - x[0];
90: PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
91: PetscCall(VecRestoreArrayRead(X, &x));
93: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
94: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
99: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
100: {
101: const PetscScalar *x;
102: PetscReal tfinal, dt;
103: User user = (User)ctx;
104: Vec interpolatedX;
106: PetscFunctionBeginUser;
107: PetscCall(TSGetTimeStep(ts, &dt));
108: PetscCall(TSGetMaxTime(ts, &tfinal));
110: while (user->next_output <= t && user->next_output <= tfinal) {
111: PetscCall(VecDuplicate(X, &interpolatedX));
112: PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
113: PetscCall(VecGetArrayRead(interpolatedX, &x));
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
115: PetscCall(VecRestoreArrayRead(interpolatedX, &x));
116: PetscCall(VecDestroy(&interpolatedX));
117: user->next_output += 0.1;
118: }
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: int main(int argc, char **argv)
123: {
124: TS ts;
125: PetscBool monitor = PETSC_FALSE;
126: PetscScalar *x_ptr;
127: PetscMPIInt size;
128: struct _n_User user;
129: PetscInt rows, cols;
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Initialize program
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: PetscFunctionBeginUser;
135: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
137: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
138: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Set runtime options
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: user.next_output = 0.0;
144: user.mu = 1.0e6;
145: user.steps = 0;
146: user.ftime = 0.5;
147: user.combined = PETSC_FALSE;
148: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
149: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
150: PetscCall(PetscOptionsGetBool(NULL, NULL, "-combined", &user.combined, NULL));
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Create necessary matrix and vectors, solve same ODE on every process
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: rows = 2;
156: cols = user.combined ? 3 : 1;
157: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jac));
158: PetscCall(MatSetSizes(user.Jac, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
159: PetscCall(MatSetFromOptions(user.Jac));
160: PetscCall(MatSetUp(user.Jac));
161: PetscCall(MatCreateVecs(user.Jac, &user.x, NULL));
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Create timestepping solver context
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
167: PetscCall(TSSetType(ts, TSBEULER));
168: PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
169: PetscCall(TSSetIJacobian(ts, user.Jac, user.Jac, IJacobian, &user));
170: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
171: PetscCall(TSSetMaxTime(ts, user.ftime));
172: if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Set initial conditions
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: PetscCall(VecGetArray(user.x, &x_ptr));
178: x_ptr[0] = 2.0;
179: x_ptr[1] = -0.66666654321;
180: PetscCall(VecRestoreArray(user.x, &x_ptr));
181: PetscCall(TSSetTimeStep(ts, 1.0 / 1024.0));
183: /* Set up forward sensitivity */
184: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
185: PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, rows, cols));
186: PetscCall(MatSetFromOptions(user.Jacp));
187: PetscCall(MatSetUp(user.Jacp));
188: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, cols, NULL, &user.sp));
189: if (user.combined) {
190: PetscCall(MatZeroEntries(user.sp));
191: PetscCall(MatShift(user.sp, 1.0));
192: } else {
193: PetscCall(MatZeroEntries(user.sp));
194: }
195: PetscCall(TSForwardSetSensitivities(ts, cols, user.sp));
196: PetscCall(TSSetRHSJacobianP(ts, user.Jacp, RHSJacobianP, &user));
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Set runtime options
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: PetscCall(TSSetFromOptions(ts));
203: PetscCall(TSSolve(ts, user.x));
204: PetscCall(TSGetSolveTime(ts, &user.ftime));
205: PetscCall(TSGetStepNumber(ts, &user.steps));
206: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, user.steps, (double)user.ftime));
207: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n ode solution \n"));
208: PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
210: if (user.combined) {
211: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n"));
212: } else {
213: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n"));
214: }
215: PetscCall(MatView(user.sp, PETSC_VIEWER_STDOUT_WORLD));
217: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: Free work space. All PETSc objects should be destroyed when they
219: are no longer needed.
220: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221: PetscCall(MatDestroy(&user.Jac));
222: PetscCall(MatDestroy(&user.sp));
223: PetscCall(MatDestroy(&user.Jacp));
224: PetscCall(VecDestroy(&user.x));
225: PetscCall(TSDestroy(&ts));
227: PetscCall(PetscFinalize());
228: return 0;
229: }
231: /*TEST
233: test:
234: args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined
235: requires: !complex !single
237: test:
238: suffix: 2
239: requires: !complex !single
240: args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5
242: TEST*/