1 #include "contexts.cxx" 2 #include "sparse.cxx" 3 #include "init.cxx" 4 #include <adolc/drivers/drivers.h> 5 #include <adolc/interfaces.h> 6 7 /* 8 REQUIRES configuration of PETSc with option --download-adolc. 9 10 For documentation on ADOL-C, see 11 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 12 */ 13 14 /* -------------------------------------------------------------------------------- 15 Drivers for RHSJacobian and IJacobian 16 ----------------------------------------------------------------------------- */ 17 18 /* 19 Compute Jacobian for explicit TS in compressed format and recover from this, using 20 precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is 21 assembled (not recommended for non-toy problems!). 22 23 Input parameters: 24 tag - tape identifier 25 u_vec - vector at which to evaluate Jacobian 26 ctx - ADOL-C context, as defined above 27 28 Output parameter: 29 A - Mat object corresponding to Jacobian 30 */ 31 PetscErrorCode PetscAdolcComputeRHSJacobian(PetscInt tag,Mat A,const PetscScalar *u_vec,void *ctx) 32 { 33 AdolcCtx *adctx = (AdolcCtx*)ctx; 34 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 35 PetscScalar **J; 36 37 PetscFunctionBegin; 38 CHKERRQ(AdolcMalloc2(m,p,&J)); 39 if (adctx->Seed) 40 fov_forward(tag,m,n,p,u_vec,adctx->Seed,NULL,J); 41 else 42 jacobian(tag,m,n,u_vec,J); 43 if (adctx->sparse) { 44 CHKERRQ(RecoverJacobian(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 45 } else { 46 for (i=0; i<m; i++) { 47 for (j=0; j<n; j++) { 48 if (fabs(J[i][j]) > 1.e-16) { 49 CHKERRQ(MatSetValues(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 50 } 51 } 52 } 53 } 54 CHKERRQ(AdolcFree2(J)); 55 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 56 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 57 PetscFunctionReturn(0); 58 } 59 60 /* 61 Compute Jacobian for explicit TS in compressed format and recover from this, using 62 precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is 63 assembled (not recommended for non-toy problems!). 64 65 Input parameters: 66 tag - tape identifier 67 u_vec - vector at which to evaluate Jacobian 68 ctx - ADOL-C context, as defined above 69 70 Output parameter: 71 A - Mat object corresponding to Jacobian 72 */ 73 PetscErrorCode PetscAdolcComputeRHSJacobianLocal(PetscInt tag,Mat A,const PetscScalar *u_vec,void *ctx) 74 { 75 AdolcCtx *adctx = (AdolcCtx*)ctx; 76 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 77 PetscScalar **J; 78 79 PetscFunctionBegin; 80 CHKERRQ(AdolcMalloc2(m,p,&J)); 81 if (adctx->Seed) 82 fov_forward(tag,m,n,p,u_vec,adctx->Seed,NULL,J); 83 else 84 jacobian(tag,m,n,u_vec,J); 85 if (adctx->sparse) { 86 CHKERRQ(RecoverJacobianLocal(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 87 } else { 88 for (i=0; i<m; i++) { 89 for (j=0; j<n; j++) { 90 if (fabs(J[i][j]) > 1.e-16) { 91 CHKERRQ(MatSetValuesLocal(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 92 } 93 } 94 } 95 } 96 CHKERRQ(AdolcFree2(J)); 97 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 98 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 99 PetscFunctionReturn(0); 100 } 101 102 /* 103 Compute Jacobian for implicit TS in compressed format and recover from this, using 104 precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is 105 assembled (not recommended for non-toy problems!). 106 107 Input parameters: 108 tag1 - tape identifier for df/dx part 109 tag2 - tape identifier for df/d(xdot) part 110 u_vec - vector at which to evaluate Jacobian 111 ctx - ADOL-C context, as defined above 112 113 Output parameter: 114 A - Mat object corresponding to Jacobian 115 */ 116 PetscErrorCode PetscAdolcComputeIJacobian(PetscInt tag1,PetscInt tag2,Mat A,const PetscScalar *u_vec,PetscReal a,void *ctx) 117 { 118 AdolcCtx *adctx = (AdolcCtx*)ctx; 119 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 120 PetscScalar **J; 121 122 PetscFunctionBegin; 123 CHKERRQ(AdolcMalloc2(m,p,&J)); 124 125 /* dF/dx part */ 126 if (adctx->Seed) 127 fov_forward(tag1,m,n,p,u_vec,adctx->Seed,NULL,J); 128 else 129 jacobian(tag1,m,n,u_vec,J); 130 CHKERRQ(MatZeroEntries(A)); 131 if (adctx->sparse) { 132 CHKERRQ(RecoverJacobian(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 133 } else { 134 for (i=0; i<m; i++) { 135 for (j=0; j<n; j++) { 136 if (fabs(J[i][j]) > 1.e-16) { 137 CHKERRQ(MatSetValues(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 138 } 139 } 140 } 141 } 142 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 143 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 144 145 /* a * dF/d(xdot) part */ 146 if (adctx->Seed) 147 fov_forward(tag2,m,n,p,u_vec,adctx->Seed,NULL,J); 148 else 149 jacobian(tag2,m,n,u_vec,J); 150 if (adctx->sparse) { 151 CHKERRQ(RecoverJacobian(A,ADD_VALUES,m,p,adctx->Rec,J,&a)); 152 } else { 153 for (i=0; i<m; i++) { 154 for (j=0; j<n; j++) { 155 if (fabs(J[i][j]) > 1.e-16) { 156 J[i][j] *= a; 157 CHKERRQ(MatSetValues(A,1,&i,1,&j,&J[i][j],ADD_VALUES)); 158 } 159 } 160 } 161 } 162 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 163 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 164 CHKERRQ(AdolcFree2(J)); 165 PetscFunctionReturn(0); 166 } 167 168 /* 169 Compute Jacobian for implicit TS in the special case where it is 170 known that the mass matrix is simply the identity. i.e. We have 171 a problem of the form 172 du/dt = F(u). 173 174 Input parameters: 175 tag - tape identifier for df/dx part 176 u_vec - vector at which to evaluate Jacobian 177 ctx - ADOL-C context, as defined above 178 179 Output parameter: 180 A - Mat object corresponding to Jacobian 181 */ 182 PetscErrorCode PetscAdolcComputeIJacobianIDMass(PetscInt tag,Mat A,PetscScalar *u_vec,PetscReal a,void *ctx) 183 { 184 AdolcCtx *adctx = (AdolcCtx*)ctx; 185 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 186 PetscScalar **J; 187 188 PetscFunctionBegin; 189 CHKERRQ(AdolcMalloc2(m,p,&J)); 190 191 /* dF/dx part */ 192 if (adctx->Seed) 193 fov_forward(tag,m,n,p,u_vec,adctx->Seed,NULL,J); 194 else 195 jacobian(tag,m,n,u_vec,J); 196 CHKERRQ(MatZeroEntries(A)); 197 if (adctx->sparse) { 198 CHKERRQ(RecoverJacobian(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 199 } else { 200 for (i=0; i<m; i++) { 201 for (j=0; j<n; j++) { 202 if (fabs(J[i][j]) > 1.e-16) { 203 CHKERRQ(MatSetValues(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 204 } 205 } 206 } 207 } 208 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 209 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 210 CHKERRQ(AdolcFree2(J)); 211 212 /* a * dF/d(xdot) part */ 213 CHKERRQ(MatShift(A,a)); 214 PetscFunctionReturn(0); 215 } 216 217 /* 218 Compute local portion of Jacobian for implicit TS in compressed format and recover from this, using 219 precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is 220 assembled (not recommended for non-toy problems!). 221 222 Input parameters: 223 tag1 - tape identifier for df/dx part 224 tag2 - tape identifier for df/d(xdot) part 225 u_vec - vector at which to evaluate Jacobian 226 ctx - ADOL-C context, as defined above 227 228 Output parameter: 229 A - Mat object corresponding to Jacobian 230 */ 231 PetscErrorCode PetscAdolcComputeIJacobianLocal(PetscInt tag1,PetscInt tag2,Mat A,PetscScalar *u_vec,PetscReal a,void *ctx) 232 { 233 AdolcCtx *adctx = (AdolcCtx*)ctx; 234 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 235 PetscScalar **J; 236 237 PetscFunctionBegin; 238 CHKERRQ(AdolcMalloc2(m,p,&J)); 239 240 /* dF/dx part */ 241 if (adctx->Seed) 242 fov_forward(tag1,m,n,p,u_vec,adctx->Seed,NULL,J); 243 else 244 jacobian(tag1,m,n,u_vec,J); 245 if (adctx->sparse) { 246 CHKERRQ(RecoverJacobianLocal(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 247 } else { 248 for (i=0; i<m; i++) { 249 for (j=0; j<n; j++) { 250 if (fabs(J[i][j]) > 1.e-16) { 251 CHKERRQ(MatSetValuesLocal(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 252 } 253 } 254 } 255 } 256 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 257 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 258 259 /* a * dF/d(xdot) part */ 260 if (adctx->Seed) 261 fov_forward(tag2,m,n,p,u_vec,adctx->Seed,NULL,J); 262 else 263 jacobian(tag2,m,n,u_vec,J); 264 if (adctx->sparse) { 265 CHKERRQ(RecoverJacobianLocal(A,ADD_VALUES,m,p,adctx->Rec,J,&a)); 266 } else { 267 for (i=0; i<m; i++) { 268 for (j=0; j<n; j++) { 269 if (fabs(J[i][j]) > 1.e-16) { 270 J[i][j] *= a; 271 CHKERRQ(MatSetValuesLocal(A,1,&i,1,&j,&J[i][j],ADD_VALUES)); 272 } 273 } 274 } 275 } 276 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 277 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 278 CHKERRQ(AdolcFree2(J)); 279 PetscFunctionReturn(0); 280 } 281 282 /* 283 Compute local portion of Jacobian for implicit TS in the special case where it is 284 known that the mass matrix is simply the identity. i.e. We have 285 a problem of the form 286 du/dt = F(u). 287 288 Input parameters: 289 tag - tape identifier for df/dx part 290 u_vec - vector at which to evaluate Jacobian 291 ctx - ADOL-C context, as defined above 292 293 Output parameter: 294 A - Mat object corresponding to Jacobian 295 */ 296 PetscErrorCode PetscAdolcComputeIJacobianLocalIDMass(PetscInt tag,Mat A,const PetscScalar *u_vec,PetscReal a,void *ctx) 297 { 298 AdolcCtx *adctx = (AdolcCtx*)ctx; 299 PetscInt i,j,m = adctx->m,n = adctx->n,p = adctx->p; 300 PetscScalar **J; 301 302 PetscFunctionBegin; 303 CHKERRQ(AdolcMalloc2(m,p,&J)); 304 305 /* dF/dx part */ 306 if (adctx->Seed) 307 fov_forward(tag,m,n,p,u_vec,adctx->Seed,NULL,J); 308 else 309 jacobian(tag,m,n,u_vec,J); 310 if (adctx->sparse) { 311 CHKERRQ(RecoverJacobianLocal(A,INSERT_VALUES,m,p,adctx->Rec,J,NULL)); 312 } else { 313 for (i=0; i<m; i++) { 314 for (j=0; j<n; j++) { 315 if (fabs(J[i][j]) > 1.e-16) { 316 CHKERRQ(MatSetValuesLocal(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 317 } 318 } 319 } 320 } 321 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 322 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 323 CHKERRQ(AdolcFree2(J)); 324 325 /* a * dF/d(xdot) part */ 326 CHKERRQ(MatShift(A,a)); 327 PetscFunctionReturn(0); 328 } 329 330 /* -------------------------------------------------------------------------------- 331 Drivers for Jacobian w.r.t. a parameter 332 ----------------------------------------------------------------------------- */ 333 334 /* 335 Compute Jacobian w.r.t a parameter for explicit TS. 336 337 Input parameters: 338 tag - tape identifier 339 u_vec - vector at which to evaluate Jacobian 340 params - the parameters w.r.t. which we differentiate 341 ctx - ADOL-C context, as defined above 342 343 Output parameter: 344 A - Mat object corresponding to Jacobian 345 */ 346 PetscErrorCode PetscAdolcComputeRHSJacobianP(PetscInt tag,Mat A,const PetscScalar *u_vec,PetscScalar *params,void *ctx) 347 { 348 AdolcCtx *adctx = (AdolcCtx*)ctx; 349 PetscInt i,j = 0,m = adctx->m,n = adctx->n,p = adctx->num_params; 350 PetscScalar **J,*concat,**S; 351 352 PetscFunctionBegin; 353 354 /* Allocate memory and concatenate independent variable values with parameter */ 355 CHKERRQ(AdolcMalloc2(m,p,&J)); 356 CHKERRQ(PetscMalloc1(n+p,&concat)); 357 CHKERRQ(AdolcMalloc2(n+p,p,&S)); 358 CHKERRQ(Subidentity(p,n,S)); 359 for (i=0; i<n; i++) concat[i] = u_vec[i]; 360 for (i=0; i<p; i++) concat[n+i] = params[i]; 361 362 /* Propagate the appropriate seed matrix through the forward mode of AD */ 363 fov_forward(tag,m,n+p,p,concat,S,NULL,J); 364 CHKERRQ(AdolcFree2(S)); 365 CHKERRQ(PetscFree(concat)); 366 367 /* Set matrix values */ 368 for (i=0; i<m; i++) { 369 for (j=0; j<p; j++) { 370 if (fabs(J[i][j]) > 1.e-16) { 371 CHKERRQ(MatSetValues(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 372 } 373 } 374 } 375 CHKERRQ(AdolcFree2(J)); 376 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 377 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 378 PetscFunctionReturn(0); 379 } 380 381 /* 382 Compute local portion of Jacobian w.r.t a parameter for explicit TS. 383 384 Input parameters: 385 tag - tape identifier 386 u_vec - vector at which to evaluate Jacobian 387 params - the parameters w.r.t. which we differentiate 388 ctx - ADOL-C context, as defined above 389 390 Output parameter: 391 A - Mat object corresponding to Jacobian 392 */ 393 PetscErrorCode PetscAdolcComputeRHSJacobianPLocal(PetscInt tag,Mat A,const PetscScalar *u_vec,PetscScalar *params,void *ctx) 394 { 395 AdolcCtx *adctx = (AdolcCtx*)ctx; 396 PetscInt i,j = 0,m = adctx->m,n = adctx->n,p = adctx->num_params; 397 PetscScalar **J,*concat,**S; 398 399 PetscFunctionBegin; 400 401 /* Allocate memory and concatenate independent variable values with parameter */ 402 CHKERRQ(AdolcMalloc2(m,p,&J)); 403 CHKERRQ(PetscMalloc1(n+p,&concat)); 404 CHKERRQ(AdolcMalloc2(n+p,p,&S)); 405 CHKERRQ(Subidentity(p,n,S)); 406 for (i=0; i<n; i++) concat[i] = u_vec[i]; 407 for (i=0; i<p; i++) concat[n+i] = params[i]; 408 409 /* Propagate the appropriate seed matrix through the forward mode of AD */ 410 fov_forward(tag,m,n+p,p,concat,S,NULL,J); 411 CHKERRQ(AdolcFree2(S)); 412 CHKERRQ(PetscFree(concat)); 413 414 /* Set matrix values */ 415 for (i=0; i<m; i++) { 416 for (j=0; j<p; j++) { 417 if (fabs(J[i][j]) > 1.e-16) { 418 CHKERRQ(MatSetValuesLocal(A,1,&i,1,&j,&J[i][j],INSERT_VALUES)); 419 } 420 } 421 } 422 CHKERRQ(AdolcFree2(J)); 423 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 424 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 425 PetscFunctionReturn(0); 426 } 427 428 /* -------------------------------------------------------------------------------- 429 Drivers for Jacobian diagonal 430 ----------------------------------------------------------------------------- */ 431 432 /* 433 Compute local portion of Jacobian diagonal for implicit TS in compressed format and recover 434 from this, using precomputed seed matrix and recovery vector. 435 436 Input parameters: 437 tag1 - tape identifier for df/dx part 438 tag2 - tape identifier for df/d(xdot) part 439 u_vec - vector at which to evaluate Jacobian 440 ctx - ADOL-C context, as defined above 441 442 Output parameter: 443 diag - Vec object corresponding to Jacobian diagonal 444 */ 445 PetscErrorCode PetscAdolcComputeIJacobianAndDiagonalLocal(PetscInt tag1,PetscInt tag2,Vec diag,PetscScalar *u_vec,PetscReal a,void *ctx) 446 { 447 AdolcCtx *adctx = (AdolcCtx*)ctx; 448 PetscInt i,m = adctx->m,n = adctx->n,p = adctx->p; 449 PetscScalar **J; 450 451 PetscFunctionBegin; 452 CHKERRQ(AdolcMalloc2(m,p,&J)); 453 454 /* dF/dx part */ 455 if (adctx->Seed) 456 fov_forward(tag1,m,n,p,u_vec,adctx->Seed,NULL,J); 457 else 458 jacobian(tag1,m,n,u_vec,J); 459 if (adctx->sparse) { 460 CHKERRQ(RecoverDiagonalLocal(diag,INSERT_VALUES,m,adctx->rec,J,NULL)); 461 } else { 462 for (i=0; i<m; i++) { 463 if (fabs(J[i][i]) > 1.e-16) { 464 CHKERRQ(VecSetValuesLocal(diag,1,&i,&J[i][i],INSERT_VALUES)); 465 } 466 } 467 } 468 CHKERRQ(VecAssemblyBegin(diag)); 469 CHKERRQ(VecAssemblyEnd(diag)); 470 471 /* a * dF/d(xdot) part */ 472 if (adctx->Seed) 473 fov_forward(tag2,m,n,p,u_vec,adctx->Seed,NULL,J); 474 else 475 jacobian(tag2,m,n,u_vec,J); 476 if (adctx->sparse) { 477 CHKERRQ(RecoverDiagonalLocal(diag,ADD_VALUES,m,adctx->rec,J,NULL)); 478 } else { 479 for (i=0; i<m; i++) { 480 if (fabs(J[i][i]) > 1.e-16) { 481 J[i][i] *= a; 482 CHKERRQ(VecSetValuesLocal(diag,1,&i,&J[i][i],ADD_VALUES)); 483 } 484 } 485 } 486 CHKERRQ(VecAssemblyBegin(diag)); 487 CHKERRQ(VecAssemblyEnd(diag)); 488 CHKERRQ(AdolcFree2(J)); 489 PetscFunctionReturn(0); 490 } 491