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] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", 113 (double)user->next_output,step,(double)t,(double)dt, 114 (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) { 173 PetscCall(TSMonitorSet(ts,Monitor,&user,NULL)); 174 } 175 176 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177 Set initial conditions 178 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179 PetscCall(VecGetArray(user.x,&x_ptr)); 180 x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321; 181 PetscCall(VecRestoreArray(user.x,&x_ptr)); 182 PetscCall(TSSetTimeStep(ts,1.0/1024.0)); 183 184 /* Set up forward sensitivity */ 185 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 186 PetscCall(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,rows,cols)); 187 PetscCall(MatSetFromOptions(user.Jacp)); 188 PetscCall(MatSetUp(user.Jacp)); 189 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,rows,cols,NULL,&user.sp)); 190 if (user.combined) { 191 PetscCall(MatZeroEntries(user.sp)); 192 PetscCall(MatShift(user.sp,1.0)); 193 } else { 194 PetscCall(MatZeroEntries(user.sp)); 195 } 196 PetscCall(TSForwardSetSensitivities(ts,cols,user.sp)); 197 PetscCall(TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user)); 198 199 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 Set runtime options 201 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202 PetscCall(TSSetFromOptions(ts)); 203 204 PetscCall(TSSolve(ts,user.x)); 205 PetscCall(TSGetSolveTime(ts,&user.ftime)); 206 PetscCall(TSGetStepNumber(ts,&user.steps)); 207 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %" PetscInt_FMT ", ftime %g\n",(double)user.mu,user.steps,(double)user.ftime)); 208 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n ode solution \n")); 209 PetscCall(VecView(user.x,PETSC_VIEWER_STDOUT_WORLD)); 210 211 if (user.combined) { 212 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n")); 213 } else { 214 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n")); 215 } 216 PetscCall(MatView(user.sp,PETSC_VIEWER_STDOUT_WORLD)); 217 218 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 219 Free work space. All PETSc objects should be destroyed when they 220 are no longer needed. 221 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 222 PetscCall(MatDestroy(&user.Jac)); 223 PetscCall(MatDestroy(&user.sp)); 224 PetscCall(MatDestroy(&user.Jacp)); 225 PetscCall(VecDestroy(&user.x)); 226 PetscCall(TSDestroy(&ts)); 227 228 PetscCall(PetscFinalize()); 229 return 0; 230 } 231 232 /*TEST 233 234 test: 235 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined 236 requires: !complex !single 237 238 test: 239 suffix: 2 240 requires: !complex !single 241 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 242 243 TEST*/ 244