1 static char help[] = "Trajectory sensitivity of a hybrid system with state-dependent switchings.\n"; 2 3 /* 4 The dynamics is described by the ODE 5 u_t = A_i u 6 7 where A_1 = [ 1 -100 8 10 1 ], 9 A_2 = [ 1 10 10 -100 1 ]. 11 The index i changes from 1 to 2 when u[1]=2.75u[0] and from 2 to 1 when u[1]=0.36u[0]. 12 Initially u=[0 1]^T and i=1. 13 14 References: 15 H. Zhang, S. Abhyankar, E. Constantinescu, M. Mihai, Discrete Adjoint Sensitivity Analysis of Hybrid Dynamical Systems With Switching, IEEE Transactions on Circuits and Systems I: Regular Papers, 64(5), May 2017 16 I. A. Hiskens, M.A. Pai, Trajectory Sensitivity Analysis of Hybrid Systems, IEEE Transactions on Circuits and Systems, Vol 47, No 2, February 2000 17 */ 18 19 #include <petscts.h> 20 21 typedef struct { 22 PetscScalar lambda1; 23 PetscScalar lambda2; 24 PetscInt mode; /* mode flag*/ 25 PetscReal print_time; 26 } AppCtx; 27 28 PetscErrorCode MyMonitor(TS ts,PetscInt stepnum,PetscReal time,Vec U,void *ctx) 29 { 30 AppCtx *actx=(AppCtx*)ctx; 31 Mat sp; 32 PetscScalar *u; 33 PetscInt nump; 34 FILE *f; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 if (time >= actx->print_time) { 39 actx->print_time += 1./256.; 40 ierr = TSForwardGetSensitivities(ts,&nump,&sp);CHKERRQ(ierr); 41 ierr = MatDenseGetColumn(sp,2,&u);CHKERRQ(ierr); 42 f = fopen("fwd_sp.out", "a"); 43 ierr = PetscFPrintf(PETSC_COMM_WORLD,f,"%20.15lf %20.15lf %20.15lf\n",time,u[0],u[1]);CHKERRQ(ierr); 44 ierr = MatDenseRestoreColumn(sp,&u);CHKERRQ(ierr); 45 fclose(f); 46 } 47 PetscFunctionReturn(0); 48 } 49 50 PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx) 51 { 52 AppCtx *actx=(AppCtx*)ctx; 53 PetscErrorCode ierr; 54 const PetscScalar *u; 55 56 PetscFunctionBegin; 57 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 58 if (actx->mode == 1) { 59 fvalue[0] = u[1]-actx->lambda1*u[0]; 60 }else if (actx->mode == 2) { 61 fvalue[0] = u[1]-actx->lambda2*u[0]; 62 } 63 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode ShiftGradients(TS ts,Vec U,AppCtx *actx) 68 { 69 Mat sp; 70 PetscScalar *x; 71 PetscScalar *u; 72 PetscErrorCode ierr; 73 PetscScalar tmp[2],A1[2][2],A2[2],denorm; 74 PetscInt nump; 75 76 PetscFunctionBegin; 77 ierr = TSForwardGetSensitivities(ts,&nump,&sp);CHKERRQ(ierr); 78 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 79 80 if (actx->mode==1) { 81 denorm = -actx->lambda1*(u[0]-100.*u[1])+1.*(10.*u[0]+u[1]); 82 A1[0][0] = 110.*u[1]*(-actx->lambda1)/denorm+1.; 83 A1[1][0] = -110.*u[0]*(-actx->lambda1)/denorm; 84 A1[0][1] = 110.*u[1]*1./denorm; 85 A1[1][1] = -110.*u[0]*1./denorm+1.; 86 87 A2[0] = 110.*u[1]*(-u[0])/denorm; 88 A2[1] = -110.*u[0]*(-u[0])/denorm; 89 }else { 90 denorm = -actx->lambda2*(u[0]+10.*u[1])+1.*(-100.*u[0]+u[1]); 91 A1[0][0] = 110.*u[1]*(actx->lambda2)/denorm+1; 92 A1[1][0] = -110.*u[0]*(actx->lambda2)/denorm; 93 A1[0][1] = -110.*u[1]*1./denorm; 94 A1[1][1] = 110.*u[0]*1./denorm+1.; 95 96 A2[0] = 0; 97 A2[1] = 0; 98 } 99 100 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 101 102 ierr = MatDenseGetColumn(sp,0,&x);CHKERRQ(ierr); 103 tmp[0] = A1[0][0]*x[0]+A1[0][1]*x[1]; 104 tmp[1] = A1[1][0]*x[0]+A1[1][1]*x[1]; 105 x[0] = tmp[0]; 106 x[1] = tmp[1]; 107 ierr = MatDenseRestoreColumn(sp,&x);CHKERRQ(ierr); 108 109 ierr = MatDenseGetColumn(sp,1,&x);CHKERRQ(ierr); 110 tmp[0] = A1[0][0]*x[0]+A1[0][1]*x[1]; 111 tmp[1] = A1[1][0]*x[0]+A1[1][1]*x[1]; 112 x[0] = tmp[0]; 113 x[1] = tmp[1]; 114 ierr = MatDenseRestoreColumn(sp,&x);CHKERRQ(ierr); 115 116 ierr = MatDenseGetColumn(sp,2,&x);CHKERRQ(ierr); 117 tmp[0] = A1[0][0]*x[0]+A1[0][1]*x[1]; 118 tmp[1] = A1[1][0]*x[0]+A1[1][1]*x[1]; 119 x[0] = tmp[0]+A2[0]; 120 x[1] = tmp[1]+A2[1]; 121 ierr = MatDenseRestoreColumn(sp,&x);CHKERRQ(ierr); 122 123 PetscFunctionReturn(0); 124 } 125 126 PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx) 127 { 128 AppCtx *actx=(AppCtx*)ctx; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 /* ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 133 ierr = ShiftGradients(ts,U,actx);CHKERRQ(ierr); 134 if (actx->mode == 1) { 135 actx->mode = 2; 136 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Change from mode 1 to 2 at t = %f \n",(double)t);CHKERRQ(ierr); */ 137 } else if (actx->mode == 2) { 138 actx->mode = 1; 139 /* ierr = PetscPrintf(PETSC_COMM_SELF,"Change from mode 2 to 1 at t = %f \n",(double)t);CHKERRQ(ierr); */ 140 } 141 PetscFunctionReturn(0); 142 } 143 144 /* 145 Defines the ODE passed to the ODE solver 146 */ 147 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 148 { 149 AppCtx *actx=(AppCtx*)ctx; 150 PetscErrorCode ierr; 151 PetscScalar *f; 152 const PetscScalar *u,*udot; 153 154 PetscFunctionBegin; 155 /* The next three lines allow us to access the entries of the vectors directly */ 156 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 157 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 158 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 159 160 if (actx->mode == 1) { 161 f[0] = udot[0]-u[0]+100*u[1]; 162 f[1] = udot[1]-10*u[0]-u[1]; 163 } else if (actx->mode == 2) { 164 f[0] = udot[0]-u[0]-10*u[1]; 165 f[1] = udot[1]+100*u[0]-u[1]; 166 } 167 168 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 169 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 170 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 /* 175 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 176 */ 177 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx) 178 { 179 AppCtx *actx=(AppCtx*)ctx; 180 PetscErrorCode ierr; 181 PetscInt rowcol[] = {0,1}; 182 PetscScalar J[2][2]; 183 const PetscScalar *u,*udot; 184 185 PetscFunctionBegin; 186 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 187 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 188 189 if (actx->mode == 1) { 190 J[0][0] = a-1; J[0][1] = 100; 191 J[1][0] = -10; J[1][1] = a-1; 192 } else if (actx->mode == 2) { 193 J[0][0] = a-1; J[0][1] = -10; 194 J[1][0] = 100; J[1][1] = a-1; 195 } 196 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 197 198 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 199 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 200 201 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 203 if (A != B) { 204 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 205 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 /* Matrix JacobianP is constant so that it only needs to be evaluated once */ 211 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat Ap,void *ctx) 212 { 213 PetscFunctionBeginUser; 214 PetscFunctionReturn(0); 215 } 216 217 int main(int argc,char **argv) 218 { 219 TS ts; /* ODE integrator */ 220 Vec U; /* solution will be stored here */ 221 Mat A; /* Jacobian matrix */ 222 Mat Ap; /* Jacobian dfdp */ 223 PetscErrorCode ierr; 224 PetscMPIInt size; 225 PetscInt n = 2; 226 PetscScalar *u; 227 AppCtx app; 228 PetscInt direction[1]; 229 PetscBool terminate[1]; 230 Mat sp; 231 PetscReal tend; 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Initialize program 234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 235 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if(ierr) return ierr; 236 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 237 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 238 app.mode = 1; 239 app.lambda1 = 2.75; 240 app.lambda2 = 0.36; 241 app.print_time = 1./256.; 242 tend = 0.125; 243 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex1fwd options","");CHKERRQ(ierr); 244 { 245 ierr = PetscOptionsReal("-lambda1","","",app.lambda1,&app.lambda1,NULL);CHKERRQ(ierr); 246 ierr = PetscOptionsReal("-lambda2","","",app.lambda2,&app.lambda2,NULL);CHKERRQ(ierr); 247 ierr = PetscOptionsReal("-tend","","",tend,&tend,NULL);CHKERRQ(ierr); 248 } 249 ierr = PetscOptionsEnd();CHKERRQ(ierr); 250 251 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 252 Create necessary matrix and vectors 253 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 254 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 255 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 256 ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr); 257 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 258 ierr = MatSetUp(A);CHKERRQ(ierr); 259 260 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 261 262 ierr = MatCreate(PETSC_COMM_WORLD,&Ap);CHKERRQ(ierr); 263 ierr = MatSetSizes(Ap,n,3,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 264 ierr = MatSetType(Ap,MATDENSE);CHKERRQ(ierr); 265 ierr = MatSetFromOptions(Ap);CHKERRQ(ierr); 266 ierr = MatSetUp(Ap);CHKERRQ(ierr); 267 ierr = MatZeroEntries(Ap);CHKERRQ(ierr); 268 269 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,3,NULL,&sp);CHKERRQ(ierr); 270 ierr = MatZeroEntries(sp);CHKERRQ(ierr); 271 ierr = MatShift(sp,1.0);CHKERRQ(ierr); 272 273 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 274 u[0] = 0; 275 u[1] = 1; 276 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 277 278 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 279 Create timestepping solver context 280 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 281 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 282 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 283 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 284 ierr = TSSetIFunction(ts,NULL,(TSIFunction)IFunction,&app);CHKERRQ(ierr); 285 ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&app);CHKERRQ(ierr); 286 287 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 288 Set initial conditions 289 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 290 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 291 292 ierr = TSForwardSetSensitivities(ts,3,sp);CHKERRQ(ierr); 293 /* Set RHS JacobianP */ 294 ierr = TSSetRHSJacobianP(ts,Ap,RHSJacobianP,&app);CHKERRQ(ierr); 295 296 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 297 Set solver options 298 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 299 ierr = TSSetMaxTime(ts,tend);CHKERRQ(ierr); 300 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 301 ierr = TSSetTimeStep(ts,1./256.);CHKERRQ(ierr); 302 ierr = TSMonitorSet(ts,MyMonitor,&app,PETSC_NULL);CHKERRQ(ierr); 303 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 304 305 /* Set directions and terminate flags for the two events */ 306 direction[0] = 0; 307 terminate[0] = PETSC_FALSE; 308 ierr = TSSetEventHandler(ts,1,direction,terminate,EventFunction,PostEventFunction,(void*)&app);CHKERRQ(ierr); 309 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 310 Run timestepping solver 311 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 312 ierr = TSSolve(ts,U);CHKERRQ(ierr); 313 314 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 315 Free work space. All PETSc objects should be destroyed when they are no longer needed. 316 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 317 ierr = MatDestroy(&A);CHKERRQ(ierr); 318 ierr = VecDestroy(&U);CHKERRQ(ierr); 319 ierr = TSDestroy(&ts);CHKERRQ(ierr); 320 321 ierr = MatDestroy(&Ap);CHKERRQ(ierr); 322 ierr = MatDestroy(&sp);CHKERRQ(ierr); 323 ierr = PetscFinalize(); 324 return(ierr); 325 } 326 327 328 /*TEST 329 330 build: 331 requires: !complex 332 333 test: 334 args: -ts_monitor 335 336 TEST*/ 337