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 User user = (User)ctx; 37 const PetscScalar *x, *xdot; 38 PetscScalar *f; 39 40 PetscFunctionBeginUser; 41 PetscCall(VecGetArrayRead(X, &x)); 42 PetscCall(VecGetArrayRead(Xdot, &xdot)); 43 PetscCall(VecGetArray(F, &f)); 44 f[0] = xdot[0] - x[1]; 45 f[1] = c21 * (xdot[0] - x[1]) + xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]); 46 PetscCall(VecRestoreArrayRead(X, &x)); 47 PetscCall(VecRestoreArrayRead(Xdot, &xdot)); 48 PetscCall(VecRestoreArray(F, &f)); 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx) { 53 User user = (User)ctx; 54 PetscInt rowcol[] = {0, 1}; 55 PetscScalar J[2][2]; 56 const PetscScalar *x; 57 58 PetscFunctionBeginUser; 59 PetscCall(VecGetArrayRead(X, &x)); 60 J[0][0] = a; 61 J[0][1] = -1.0; 62 J[1][0] = c21 * a + user->mu * (1.0 + 2.0 * x[0] * x[1]); 63 J[1][1] = -c21 + a - user->mu * (1.0 - x[0] * x[0]); 64 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 65 PetscCall(VecRestoreArrayRead(X, &x)); 66 67 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 68 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 69 if (A != B) { 70 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 71 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 72 } 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx) { 77 User user = (User)ctx; 78 PetscInt row[] = {0, 1}, col[] = {0}; 79 PetscScalar J[2][1]; 80 const PetscScalar *x; 81 82 PetscFunctionBeginUser; 83 if (user->combined) col[0] = 2; 84 PetscCall(VecGetArrayRead(X, &x)); 85 J[0][0] = 0; 86 J[1][0] = (1. - x[0] * x[0]) * x[1] - x[0]; 87 PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES)); 88 PetscCall(VecRestoreArrayRead(X, &x)); 89 90 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 91 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 92 PetscFunctionReturn(0); 93 } 94 95 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 96 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx) { 97 const PetscScalar *x; 98 PetscReal tfinal, dt; 99 User user = (User)ctx; 100 Vec interpolatedX; 101 102 PetscFunctionBeginUser; 103 PetscCall(TSGetTimeStep(ts, &dt)); 104 PetscCall(TSGetMaxTime(ts, &tfinal)); 105 106 while (user->next_output <= t && user->next_output <= tfinal) { 107 PetscCall(VecDuplicate(X, &interpolatedX)); 108 PetscCall(TSInterpolate(ts, user->next_output, interpolatedX)); 109 PetscCall(VecGetArrayRead(interpolatedX, &x)); 110 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]))); 111 PetscCall(VecRestoreArrayRead(interpolatedX, &x)); 112 PetscCall(VecDestroy(&interpolatedX)); 113 user->next_output += 0.1; 114 } 115 PetscFunctionReturn(0); 116 } 117 118 int main(int argc, char **argv) { 119 TS ts; 120 PetscBool monitor = PETSC_FALSE; 121 PetscScalar *x_ptr; 122 PetscMPIInt size; 123 struct _n_User user; 124 PetscInt rows, cols; 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Initialize program 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscFunctionBeginUser; 130 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 131 132 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 133 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Set runtime options 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 user.next_output = 0.0; 139 user.mu = 1.0e6; 140 user.steps = 0; 141 user.ftime = 0.5; 142 user.combined = PETSC_FALSE; 143 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL)); 144 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL)); 145 PetscCall(PetscOptionsGetBool(NULL, NULL, "-combined", &user.combined, NULL)); 146 147 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 148 Create necessary matrix and vectors, solve same ODE on every process 149 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 150 rows = 2; 151 cols = user.combined ? 3 : 1; 152 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jac)); 153 PetscCall(MatSetSizes(user.Jac, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 154 PetscCall(MatSetFromOptions(user.Jac)); 155 PetscCall(MatSetUp(user.Jac)); 156 PetscCall(MatCreateVecs(user.Jac, &user.x, NULL)); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Create timestepping solver context 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 162 PetscCall(TSSetType(ts, TSBEULER)); 163 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 164 PetscCall(TSSetIJacobian(ts, user.Jac, user.Jac, IJacobian, &user)); 165 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 166 PetscCall(TSSetMaxTime(ts, user.ftime)); 167 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL)); 168 169 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 Set initial conditions 171 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 172 PetscCall(VecGetArray(user.x, &x_ptr)); 173 x_ptr[0] = 2.0; 174 x_ptr[1] = -0.66666654321; 175 PetscCall(VecRestoreArray(user.x, &x_ptr)); 176 PetscCall(TSSetTimeStep(ts, 1.0 / 1024.0)); 177 178 /* Set up forward sensitivity */ 179 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp)); 180 PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, rows, cols)); 181 PetscCall(MatSetFromOptions(user.Jacp)); 182 PetscCall(MatSetUp(user.Jacp)); 183 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, cols, NULL, &user.sp)); 184 if (user.combined) { 185 PetscCall(MatZeroEntries(user.sp)); 186 PetscCall(MatShift(user.sp, 1.0)); 187 } else { 188 PetscCall(MatZeroEntries(user.sp)); 189 } 190 PetscCall(TSForwardSetSensitivities(ts, cols, user.sp)); 191 PetscCall(TSSetRHSJacobianP(ts, user.Jacp, RHSJacobianP, &user)); 192 193 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 Set runtime options 195 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196 PetscCall(TSSetFromOptions(ts)); 197 198 PetscCall(TSSolve(ts, user.x)); 199 PetscCall(TSGetSolveTime(ts, &user.ftime)); 200 PetscCall(TSGetStepNumber(ts, &user.steps)); 201 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, user.steps, (double)user.ftime)); 202 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n ode solution \n")); 203 PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 204 205 if (user.combined) { 206 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n")); 207 } else { 208 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n")); 209 } 210 PetscCall(MatView(user.sp, PETSC_VIEWER_STDOUT_WORLD)); 211 212 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213 Free work space. All PETSc objects should be destroyed when they 214 are no longer needed. 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 PetscCall(MatDestroy(&user.Jac)); 217 PetscCall(MatDestroy(&user.sp)); 218 PetscCall(MatDestroy(&user.Jacp)); 219 PetscCall(VecDestroy(&user.x)); 220 PetscCall(TSDestroy(&ts)); 221 222 PetscCall(PetscFinalize()); 223 return 0; 224 } 225 226 /*TEST 227 228 test: 229 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined 230 requires: !complex !single 231 232 test: 233 suffix: 2 234 requires: !complex !single 235 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 236 237 TEST*/ 238