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; J[0][1] = -1.0; 63 J[1][0] = c21*a + user->mu*(1.0 + 2.0*x[0]*x[1]); 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 { 78 User user = (User)ctx; 79 PetscInt row[] = {0,1},col[]={0}; 80 PetscScalar J[2][1]; 81 const PetscScalar *x; 82 83 PetscFunctionBeginUser; 84 if (user->combined) col[0] = 2; 85 PetscCall(VecGetArrayRead(X,&x)); 86 J[0][0] = 0; 87 J[1][0] = (1.-x[0]*x[0])*x[1]-x[0]; 88 PetscCall(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 89 PetscCall(VecRestoreArrayRead(X,&x)); 90 91 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 93 PetscFunctionReturn(0); 94 } 95 96 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 97 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx) 98 { 99 const PetscScalar *x; 100 PetscReal tfinal, dt; 101 User user = (User)ctx; 102 Vec interpolatedX; 103 104 PetscFunctionBeginUser; 105 PetscCall(TSGetTimeStep(ts,&dt)); 106 PetscCall(TSGetMaxTime(ts,&tfinal)); 107 108 while (user->next_output <= t && user->next_output <= tfinal) { 109 PetscCall(VecDuplicate(X,&interpolatedX)); 110 PetscCall(TSInterpolate(ts,user->next_output,interpolatedX)); 111 PetscCall(VecGetArrayRead(interpolatedX,&x)); 112 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", 113 user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]), 114 (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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 135 136 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 137 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 138 139 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 140 Set runtime options 141 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 142 user.next_output = 0.0; 143 user.mu = 1.0e6; 144 user.steps = 0; 145 user.ftime = 0.5; 146 user.combined = PETSC_FALSE; 147 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 148 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 149 PetscCall(PetscOptionsGetBool(NULL,NULL,"-combined",&user.combined,NULL)); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Create necessary matrix and vectors, solve same ODE on every process 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 rows = 2; 155 cols = user.combined ? 3 : 1; 156 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Jac)); 157 PetscCall(MatSetSizes(user.Jac,PETSC_DECIDE,PETSC_DECIDE,2,2)); 158 PetscCall(MatSetFromOptions(user.Jac)); 159 PetscCall(MatSetUp(user.Jac)); 160 PetscCall(MatCreateVecs(user.Jac,&user.x,NULL)); 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Create timestepping solver context 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 166 PetscCall(TSSetType(ts,TSBEULER)); 167 PetscCall(TSSetIFunction(ts,NULL,IFunction,&user)); 168 PetscCall(TSSetIJacobian(ts,user.Jac,user.Jac,IJacobian,&user)); 169 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 170 PetscCall(TSSetMaxTime(ts,user.ftime)); 171 if (monitor) { 172 PetscCall(TSMonitorSet(ts,Monitor,&user,NULL)); 173 } 174 175 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 176 Set initial conditions 177 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 178 PetscCall(VecGetArray(user.x,&x_ptr)); 179 x_ptr[0] = 2.0; 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 %D, 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