1 static char help[] = "Performs adjoint sensitivity analysis for a van der Pol like \n\ 2 equation with time dependent parameters using two approaches : \n\ 3 track : track only local sensitivities at each adjoint step \n\ 4 and accumulate them in a global array \n\ 5 global : track parameters at all timesteps together \n\ 6 Choose one of the two at runtime by -sa_method {track,global}. \n"; 7 8 /* 9 Concepts: TS^adjoint for time dependent parameters 10 Concepts: TS^Customized adjoint monitor based sensitivity tracking 11 Concepts: TS^All at once approach to sensitivity tracking 12 Processors: 1 13 */ 14 15 /* 16 Simple example to demonstrate TSAdjoint capabilities for time dependent params 17 without integral cost terms using either a tracking or global method. 18 19 Modify the Van Der Pol Eq to : 20 [u1'] = [mu1(t)*u1] 21 [u2'] = [mu2(t)*((1-u1^2)*u2-u1)] 22 (with initial conditions & params independent) 23 24 Define uref to be solution with initail conditions (2,-2/3), mu=(1,1e3) 25 - u_ref : (1.5967,-1.02969) 26 27 Define const function as cost = 2-norm(u - u_ref); 28 29 Initialization for the adjoint TS : 30 - dcost/dy|final_time = 2*(u-u_ref)|final_time 31 - dcost/dp|final_time = 0 32 33 The tracking method only tracks local sensitivity at each time step 34 and accumulates these sensitivities in a global array. Since the structure 35 of the equations being solved at each time step does not change, the jacobian 36 wrt parameters is defined analogous to constant RHSJacobian for a liner 37 TSSolve and the size of the jacP is independent of the number of time 38 steps. Enable this mode of adjoint analysis by -sa_method track. 39 40 The global method combines the parameters at all timesteps and tracks them 41 together. Thus, the columns of the jacP matrix are filled dependent upon the 42 time step. Also, the dimensions of the jacP matrix now depend upon the number 43 of time steps. Enable this mode of adjoint analysis by -sa_method global. 44 45 Since the equations here have parameters at predefined time steps, this 46 example should be run with non adaptive time stepping solvers only. This 47 can be ensured by -ts_adapt_type none (which is the default behavior only 48 for certain TS solvers like TSCN. If using an explicit TS like TSRK, 49 please be sure to add the aforementioned option to disable adaptive 50 timestepping.) 51 */ 52 53 /* 54 Include "petscts.h" so that we can use TS solvers. Note that this file 55 automatically includes: 56 petscsys.h - base PETSc routines petscvec.h - vectors 57 petscmat.h - matrices 58 petscis.h - index sets petscksp.h - Krylov subspace methods 59 petscviewer.h - viewers petscpc.h - preconditioners 60 petscksp.h - linear solvers petscsnes.h - nonlinear solvers 61 */ 62 #include <petscts.h> 63 64 extern PetscErrorCode RHSFunction(TS ,PetscReal ,Vec ,Vec ,void *); 65 extern PetscErrorCode RHSJacobian(TS ,PetscReal ,Vec ,Mat ,Mat ,void *); 66 extern PetscErrorCode RHSJacobianP_track(TS ,PetscReal ,Vec ,Mat ,void *); 67 extern PetscErrorCode RHSJacobianP_global(TS ,PetscReal ,Vec ,Mat ,void *); 68 extern PetscErrorCode Monitor(TS ,PetscInt ,PetscReal ,Vec ,void *); 69 extern PetscErrorCode AdjointMonitor(TS ,PetscInt ,PetscReal ,Vec ,PetscInt ,Vec *, Vec *,void *); 70 71 /* 72 User-defined application context - contains data needed by the 73 application-provided call-back routines. 74 */ 75 76 typedef struct { 77 /*------------- Forward solve data structures --------------*/ 78 PetscInt max_steps; /* number of steps to run ts for */ 79 PetscReal final_time; /* final time to integrate to*/ 80 PetscReal time_step; /* ts integration time step */ 81 Vec mu1; /* time dependent params */ 82 Vec mu2; /* time dependent params */ 83 Vec U; /* solution vector */ 84 Mat A; /* Jacobian matrix */ 85 86 /*------------- Adjoint solve data structures --------------*/ 87 Mat Jacp; /* JacobianP matrix */ 88 Vec lambda; /* adjoint variable */ 89 Vec mup; /* adjoint variable */ 90 91 /*------------- Global accumation vecs for monitor based tracking --------------*/ 92 Vec sens_mu1; /* global sensitivity accumulation */ 93 Vec sens_mu2; /* global sensitivity accumulation */ 94 PetscInt adj_idx; /* to keep track of adjoint solve index */ 95 } AppCtx; 96 97 typedef enum {SA_TRACK, SA_GLOBAL} SAMethod; 98 static const char *const SAMethods[] = {"TRACK","GLOBAL","SAMethod","SA_",0}; 99 100 /* ----------------------- Explicit form of the ODE -------------------- */ 101 102 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 103 { 104 PetscErrorCode ierr; 105 AppCtx *user = (AppCtx*) ctx; 106 PetscScalar *f; 107 PetscInt curr_step; 108 const PetscScalar *u; 109 const PetscScalar *mu1; 110 const PetscScalar *mu2; 111 112 PetscFunctionBeginUser; 113 ierr = TSGetStepNumber(ts,&curr_step);CHKERRQ(ierr); 114 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 115 ierr = VecGetArrayRead(user->mu1,&mu1);CHKERRQ(ierr); 116 ierr = VecGetArrayRead(user->mu2,&mu2);CHKERRQ(ierr); 117 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 118 f[0] = mu1[curr_step]*u[1]; 119 f[1] = mu2[curr_step]*((1.-u[0]*u[0])*u[1]-u[0]); 120 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 121 ierr = VecRestoreArrayRead(user->mu1,&mu1);CHKERRQ(ierr); 122 ierr = VecRestoreArrayRead(user->mu2,&mu2);CHKERRQ(ierr); 123 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 128 { 129 PetscErrorCode ierr; 130 AppCtx *user = (AppCtx*) ctx; 131 PetscInt rowcol[] = {0,1}; 132 PetscScalar J[2][2]; 133 PetscInt curr_step; 134 const PetscScalar *u; 135 const PetscScalar *mu1; 136 const PetscScalar *mu2; 137 138 PetscFunctionBeginUser; 139 ierr = TSGetStepNumber(ts,&curr_step);CHKERRQ(ierr); 140 ierr = VecGetArrayRead(user->mu1,&mu1);CHKERRQ(ierr); 141 ierr = VecGetArrayRead(user->mu2,&mu2);CHKERRQ(ierr); 142 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 143 J[0][0] = 0; 144 J[1][0] = -mu2[curr_step]*(2.0*u[1]*u[0]+1.); 145 J[0][1] = mu1[curr_step]; 146 J[1][1] = mu2[curr_step]*(1.0-u[0]*u[0]); 147 ierr = MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 148 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 149 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 151 ierr = VecRestoreArrayRead(user->mu1,&mu1);CHKERRQ(ierr); 152 ierr = VecRestoreArrayRead(user->mu2,&mu2);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 /* ------------------ Jacobian wrt parameters for tracking method ------------------ */ 157 158 PetscErrorCode RHSJacobianP_track(TS ts,PetscReal t,Vec U,Mat A,void *ctx) 159 { 160 PetscErrorCode ierr; 161 PetscInt row[] = {0,1},col[] = {0,1}; 162 PetscScalar J[2][2]; 163 const PetscScalar *u; 164 165 PetscFunctionBeginUser; 166 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 167 J[0][0] = u[1]; 168 J[1][0] = 0; 169 J[0][1] = 0; 170 J[1][1] = (1.-u[0]*u[0])*u[1]-u[0]; 171 ierr = MatSetValues(A,2,row,2,col,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 172 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 173 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 174 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 178 /* ------------------ Jacobian wrt parameters for global method ------------------ */ 179 180 PetscErrorCode RHSJacobianP_global(TS ts,PetscReal t,Vec U,Mat A,void *ctx) 181 { 182 PetscErrorCode ierr; 183 PetscInt row[] = {0,1},col[] = {0,1}; 184 PetscScalar J[2][2]; 185 const PetscScalar *u; 186 PetscInt curr_step; 187 188 PetscFunctionBeginUser; 189 ierr = TSGetStepNumber(ts,&curr_step);CHKERRQ(ierr); 190 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 191 J[0][0] = u[1]; 192 J[1][0] = 0; 193 J[0][1] = 0; 194 J[1][1] = (1.-u[0]*u[0])*u[1]-u[0]; 195 col[0] = (curr_step)*2; 196 col[1] = (curr_step)*2+1; 197 ierr = MatSetValues(A,2,row,2,col,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 198 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 199 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 200 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 /* Dump solution to console if called */ 205 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx) 206 { 207 PetscErrorCode ierr; 208 209 PetscFunctionBeginUser; 210 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n Solution at time %e is \n", t);CHKERRQ(ierr); 211 ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214 215 /* Customized adjoint monitor to keep track of local 216 sensitivities by storing them in a global sensitivity array. 217 Note : This routine is only used for the tracking method. */ 218 PetscErrorCode AdjointMonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *ctx) 219 { 220 PetscErrorCode ierr; 221 AppCtx *user = (AppCtx*) ctx; 222 PetscInt curr_step; 223 PetscScalar *sensmu1_glob; 224 PetscScalar *sensmu2_glob; 225 const PetscScalar *sensmu_loc; 226 227 PetscFunctionBeginUser; 228 ierr = TSGetStepNumber(ts,&curr_step);CHKERRQ(ierr); 229 /* Note that we skip the first call to the monitor in the adjoint 230 solve since the sensitivities are already set (during 231 initialization of adjoint vectors). 232 We also note that each indvidial TSAdjointSolve calls the monitor 233 twice, once at the step it is integrating from and once at the step 234 it integrates to. Only the second call is useful for transferring 235 local sensitivities to the global array. */ 236 if (curr_step == user->adj_idx) { 237 PetscFunctionReturn(0); 238 } else { 239 ierr = VecGetArrayRead(*mu,&sensmu_loc);CHKERRQ(ierr); 240 ierr = VecGetArray(user->sens_mu1,&sensmu1_glob);CHKERRQ(ierr); 241 ierr = VecGetArray(user->sens_mu2,&sensmu2_glob);CHKERRQ(ierr); 242 sensmu1_glob[curr_step] = sensmu_loc[0]; 243 sensmu2_glob[curr_step] = sensmu_loc[1]; 244 ierr = VecRestoreArray(user->sens_mu1,&sensmu1_glob);CHKERRQ(ierr); 245 ierr = VecRestoreArray(user->sens_mu2,&sensmu2_glob);CHKERRQ(ierr); 246 ierr = VecRestoreArrayRead(*mu,&sensmu_loc);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 } 250 251 int main(int argc,char **argv) 252 { 253 TS ts; 254 AppCtx user; 255 PetscScalar *x_ptr,*y_ptr,*u_ptr; 256 PetscMPIInt size; 257 PetscBool monitor = PETSC_FALSE; 258 SAMethod sa = SA_GLOBAL; 259 PetscErrorCode ierr; 260 261 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 262 Initialize program 263 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 264 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 265 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 266 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 267 268 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 269 Set runtime options 270 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 271 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"SA analysis options.","");CHKERRQ(ierr);{ 272 ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 273 ierr = PetscOptionsEnum("-sa_method","Sensitivity analysis method (track or global)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL);CHKERRQ(ierr); 274 } 275 ierr = PetscOptionsEnd();CHKERRQ(ierr); 276 277 user.final_time = 0.1; 278 user.max_steps = 5; 279 user.time_step = user.final_time/user.max_steps; 280 281 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 282 Create necessary matrix and vectors for forward solve. 283 Create Jacp matrix for adjoint solve. 284 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 285 ierr = VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.mu1);CHKERRQ(ierr); 286 ierr = VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.mu2);CHKERRQ(ierr); 287 ierr = VecSet(user.mu1,1.25);CHKERRQ(ierr); 288 ierr = VecSet(user.mu2,1.0e2);CHKERRQ(ierr); 289 290 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 291 For tracking method : create the global sensitivity array to 292 accumulate sensitivity with respect to parameters at each step. 293 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 294 if (sa == SA_TRACK) { 295 ierr = VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.sens_mu1);CHKERRQ(ierr); 296 ierr = VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.sens_mu2);CHKERRQ(ierr); 297 } 298 299 ierr = MatCreate(PETSC_COMM_WORLD,&user.A);CHKERRQ(ierr); 300 ierr = MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 301 ierr = MatSetFromOptions(user.A);CHKERRQ(ierr); 302 ierr = MatSetUp(user.A);CHKERRQ(ierr); 303 ierr = MatCreateVecs(user.A,&user.U,NULL);CHKERRQ(ierr); 304 305 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 306 Note that the dimensions of the Jacp matrix depend upon the 307 sensitivity analysis method being used ! 308 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 309 ierr = MatCreate(PETSC_COMM_WORLD,&user.Jacp);CHKERRQ(ierr); 310 if (sa == SA_TRACK) { 311 ierr = MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 312 } 313 if (sa == SA_GLOBAL) { 314 ierr = MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,user.max_steps*2);CHKERRQ(ierr); 315 } 316 ierr = MatSetFromOptions(user.Jacp);CHKERRQ(ierr); 317 ierr = MatSetUp(user.Jacp);CHKERRQ(ierr); 318 319 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 320 Create timestepping solver context 321 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 322 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 323 ierr = TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); 324 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 325 326 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 327 ierr = TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user);CHKERRQ(ierr); 328 if (sa == SA_TRACK) { 329 ierr = TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP_track,&user);CHKERRQ(ierr); 330 } 331 if (sa == SA_GLOBAL) { 332 ierr = TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP_global,&user);CHKERRQ(ierr); 333 } 334 335 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 336 ierr = TSSetMaxTime(ts,user.final_time);CHKERRQ(ierr); 337 ierr = TSSetTimeStep(ts,user.final_time/user.max_steps);CHKERRQ(ierr); 338 339 if (monitor) { 340 ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr); 341 } 342 if (sa == SA_TRACK) { 343 ierr = TSAdjointMonitorSet(ts,AdjointMonitor,&user,NULL);CHKERRQ(ierr); 344 } 345 346 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 347 Set initial conditions 348 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 349 ierr = VecGetArray(user.U,&x_ptr);CHKERRQ(ierr); 350 x_ptr[0] = 2.0; 351 x_ptr[1] = -2.0/3.0; 352 ierr = VecRestoreArray(user.U,&x_ptr);CHKERRQ(ierr); 353 354 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 355 Save trajectory of solution so that TSAdjointSolve() may be used 356 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 357 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr); 358 359 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 360 Set runtime options 361 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 362 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 363 364 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 365 Execute forward model and print solution. 366 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 367 ierr = TSSolve(ts,user.U);CHKERRQ(ierr); 368 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n Solution of forward TS :\n");CHKERRQ(ierr); 369 ierr = VecView(user.U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 370 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n Forward TS solve successfull! Adjoint run begins!\n");CHKERRQ(ierr); 371 372 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 373 Adjoint model starts here! Create adjoint vectors. 374 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 375 ierr = MatCreateVecs(user.A,&user.lambda,NULL);CHKERRQ(ierr); 376 ierr = MatCreateVecs(user.Jacp,&user.mup,NULL);CHKERRQ(ierr); 377 378 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 379 Set initial conditions for the adjoint vector 380 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 381 ierr = VecGetArray(user.U,&u_ptr);CHKERRQ(ierr); 382 ierr = VecGetArray(user.lambda,&y_ptr);CHKERRQ(ierr); 383 y_ptr[0] = 2*(u_ptr[0] - 1.5967); 384 y_ptr[1] = 2*(u_ptr[1] - -(1.02969)); 385 ierr = VecRestoreArray(user.lambda,&y_ptr);CHKERRQ(ierr); 386 ierr = VecRestoreArray(user.U,&y_ptr);CHKERRQ(ierr); 387 ierr = VecSet(user.mup,0);CHKERRQ(ierr); 388 389 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 390 Set number of cost functions. 391 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 392 ierr = TSSetCostGradients(ts,1,&user.lambda,&user.mup);CHKERRQ(ierr); 393 394 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 395 The adjoint vector mup has to be reset for each adjoint step when 396 using the tracking method as we want to treat the parameters at each 397 time step one at a time and prevent accumulation of the sensitivities 398 from parameters at previous time steps. 399 This is not necessary for the global method as each time dependent 400 parameter is treated as an independent parameter. 401 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 402 if (sa == SA_TRACK) { 403 for (user.adj_idx=user.max_steps; user.adj_idx>0; user.adj_idx--) { 404 ierr = VecSet(user.mup,0);CHKERRQ(ierr); 405 ierr = TSAdjointSetSteps(ts, 1);CHKERRQ(ierr); 406 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 407 } 408 } 409 if (sa == SA_GLOBAL) { 410 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 411 } 412 413 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 414 Dispaly adjoint sensitivities wrt parameters and initial conditions 415 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 416 if (sa == SA_TRACK) { 417 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt mu1: d[cost]/d[mu1]\n");CHKERRQ(ierr); 418 ierr = VecView(user.sens_mu1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 419 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt mu2: d[cost]/d[mu2]\n");CHKERRQ(ierr); 420 ierr = VecView(user.sens_mu2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 421 } 422 423 if (sa == SA_GLOBAL) { 424 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt params: d[cost]/d[p], where p refers to \n\ 425 the interlaced vector made by combining mu1,mu2\n");CHKERRQ(ierr); 426 ierr = VecView(user.mup,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 427 } 428 429 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[cost]/d[u(t=0)]\n");CHKERRQ(ierr); 430 ierr = VecView(user.lambda,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 431 432 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 433 Free work space! 434 All PETSc objects should be destroyed when they are no longer needed. 435 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 436 ierr = MatDestroy(&user.A);CHKERRQ(ierr); 437 ierr = MatDestroy(&user.Jacp);CHKERRQ(ierr); 438 ierr = VecDestroy(&user.U);CHKERRQ(ierr); 439 ierr = VecDestroy(&user.lambda);CHKERRQ(ierr); 440 ierr = VecDestroy(&user.mup);CHKERRQ(ierr); 441 ierr = VecDestroy(&user.mu1);CHKERRQ(ierr); 442 ierr = VecDestroy(&user.mu2);CHKERRQ(ierr); 443 if (sa == SA_TRACK) { 444 ierr = VecDestroy(&user.sens_mu1);CHKERRQ(ierr); 445 ierr = VecDestroy(&user.sens_mu2);CHKERRQ(ierr); 446 } 447 ierr = TSDestroy(&ts);CHKERRQ(ierr); 448 ierr = PetscFinalize(); 449 return(ierr); 450 } 451 452 /*TEST 453 454 test: 455 requires: !complex 456 suffix : track 457 args : -sa_method track 458 459 test: 460 requires: !complex 461 suffix : global 462 args : -sa_method global 463 464 TEST*/ 465 466