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