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"; 7 8 /* ------------------------------------------------------------------------ 9 10 This code demonstrates how to compute trajectory sensitivties w.r.t. the stiffness parameter mu. 11 1) Use two vectors s and sp for sensitivities w.r.t. initial values and paraeters 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' 13 14 ------------------------------------------------------------------------- */ 15 #include <petscts.h> 16 #include <petsctao.h> 17 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 }; 31 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; 40 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(0); 51 } 52 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; 59 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)); 68 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(0); 76 } 77 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; 84 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)); 92 93 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 94 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 95 PetscFunctionReturn(0); 96 } 97 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; 105 106 PetscFunctionBeginUser; 107 PetscCall(TSGetTimeStep(ts, &dt)); 108 PetscCall(TSGetMaxTime(ts, &tfinal)); 109 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(0); 120 } 121 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; 130 131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 132 Initialize program 133 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 134 PetscFunctionBeginUser; 135 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 136 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!"); 139 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)); 151 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)); 162 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)); 173 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)); 182 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)); 197 198 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 199 Set runtime options 200 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 201 PetscCall(TSSetFromOptions(ts)); 202 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)); 209 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)); 216 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)); 226 227 PetscCall(PetscFinalize()); 228 return 0; 229 } 230 231 /*TEST 232 233 test: 234 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined 235 requires: !complex !single 236 237 test: 238 suffix: 2 239 requires: !complex !single 240 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 241 242 TEST*/ 243