1 #include <petsc/private/dmimpl.h> 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 4 typedef struct { 5 PetscErrorCode (*boundarylocal)(DM, PetscReal, Vec, Vec, void *); 6 PetscErrorCode (*ifunctionlocal)(DM, PetscReal, Vec, Vec, Vec, void *); 7 PetscErrorCode (*ijacobianlocal)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 8 PetscErrorCode (*rhsfunctionlocal)(DM, PetscReal, Vec, Vec, void *); 9 void *boundarylocalctx; 10 void *ifunctionlocalctx; 11 void *ijacobianlocalctx; 12 void *rhsfunctionlocalctx; 13 Vec lumpedmassinv; 14 Mat mass; 15 KSP kspmass; 16 } DMTS_Local; 17 18 static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm) 19 { 20 PetscFunctionBegin; 21 PetscCall(PetscFree(tdm->data)); 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm) 26 { 27 PetscFunctionBegin; 28 PetscCall(PetscNew((DMTS_Local **)&tdm->data)); 29 if (oldtdm->data) PetscCall(PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local))); 30 PetscFunctionReturn(PETSC_SUCCESS); 31 } 32 33 static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts) 34 { 35 PetscFunctionBegin; 36 *dmlocalts = NULL; 37 if (!tdm->data) { 38 PetscCall(PetscNew((DMTS_Local **)&tdm->data)); 39 40 tdm->ops->destroy = DMTSDestroy_DMLocal; 41 tdm->ops->duplicate = DMTSDuplicate_DMLocal; 42 } 43 *dmlocalts = (DMTS_Local *)tdm->data; 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx) 48 { 49 DM dm; 50 Vec locX, locX_t, locF; 51 DMTS_Local *dmlocalts = (DMTS_Local *)ctx; 52 53 PetscFunctionBegin; 54 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 55 PetscValidHeaderSpecific(X, VEC_CLASSID, 3); 56 PetscValidHeaderSpecific(X_t, VEC_CLASSID, 4); 57 PetscValidHeaderSpecific(F, VEC_CLASSID, 5); 58 PetscCall(TSGetDM(ts, &dm)); 59 PetscCall(DMGetLocalVector(dm, &locX)); 60 PetscCall(DMGetLocalVector(dm, &locX_t)); 61 PetscCall(DMGetLocalVector(dm, &locF)); 62 PetscCall(VecZeroEntries(locX)); 63 PetscCall(VecZeroEntries(locX_t)); 64 if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx)); 65 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX)); 66 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX)); 67 PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t)); 68 PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t)); 69 PetscCall(VecZeroEntries(locF)); 70 CHKMEMQ; 71 PetscCall((*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx)); 72 CHKMEMQ; 73 PetscCall(VecZeroEntries(F)); 74 PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F)); 75 PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F)); 76 PetscCall(DMRestoreLocalVector(dm, &locX)); 77 PetscCall(DMRestoreLocalVector(dm, &locX_t)); 78 PetscCall(DMRestoreLocalVector(dm, &locF)); 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx) 83 { 84 DM dm; 85 Vec locX, locF; 86 DMTS_Local *dmlocalts = (DMTS_Local *)ctx; 87 88 PetscFunctionBegin; 89 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 90 PetscValidHeaderSpecific(X, VEC_CLASSID, 3); 91 PetscValidHeaderSpecific(F, VEC_CLASSID, 4); 92 PetscCall(TSGetDM(ts, &dm)); 93 PetscCall(DMGetLocalVector(dm, &locX)); 94 PetscCall(DMGetLocalVector(dm, &locF)); 95 PetscCall(VecZeroEntries(locX)); 96 if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx)); 97 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX)); 98 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX)); 99 PetscCall(VecZeroEntries(locF)); 100 CHKMEMQ; 101 PetscCall((*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx)); 102 CHKMEMQ; 103 PetscCall(VecZeroEntries(F)); 104 PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F)); 105 PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F)); 106 if (dmlocalts->lumpedmassinv) { 107 PetscCall(VecPointwiseMult(F, dmlocalts->lumpedmassinv, F)); 108 } else if (dmlocalts->kspmass) { 109 Vec tmp; 110 111 PetscCall(DMGetGlobalVector(dm, &tmp)); 112 PetscCall(KSPSolve(dmlocalts->kspmass, F, tmp)); 113 PetscCall(VecCopy(tmp, F)); 114 PetscCall(DMRestoreGlobalVector(dm, &tmp)); 115 } 116 PetscCall(DMRestoreLocalVector(dm, &locX)); 117 PetscCall(DMRestoreLocalVector(dm, &locF)); 118 PetscFunctionReturn(PETSC_SUCCESS); 119 } 120 121 static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx) 122 { 123 DM dm; 124 Vec locX, locX_t; 125 DMTS_Local *dmlocalts = (DMTS_Local *)ctx; 126 127 PetscFunctionBegin; 128 PetscCall(TSGetDM(ts, &dm)); 129 if (dmlocalts->ijacobianlocal) { 130 PetscCall(DMGetLocalVector(dm, &locX)); 131 PetscCall(DMGetLocalVector(dm, &locX_t)); 132 PetscCall(VecZeroEntries(locX)); 133 PetscCall(VecZeroEntries(locX_t)); 134 if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx)); 135 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX)); 136 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX)); 137 PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t)); 138 PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t)); 139 CHKMEMQ; 140 PetscCall((*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx)); 141 CHKMEMQ; 142 PetscCall(DMRestoreLocalVector(dm, &locX)); 143 PetscCall(DMRestoreLocalVector(dm, &locX_t)); 144 } else { 145 MatFDColoring fdcoloring; 146 PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring)); 147 if (!fdcoloring) { 148 ISColoring coloring; 149 150 PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring)); 151 PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring)); 152 PetscCall(ISColoringDestroy(&coloring)); 153 switch (dm->coloringtype) { 154 case IS_COLORING_GLOBAL: 155 PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))TSComputeIFunction_DMLocal, dmlocalts)); 156 break; 157 default: 158 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]); 159 } 160 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix)); 161 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 162 PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring)); 163 PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring)); 164 PetscCall(PetscObjectDereference((PetscObject)fdcoloring)); 165 166 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 167 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 168 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 169 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 170 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 171 */ 172 PetscCall(PetscObjectDereference((PetscObject)dm)); 173 } 174 PetscCall(MatFDColoringApply(B, fdcoloring, X, ts)); 175 } 176 /* This will be redundant if the user called both, but it's too common to forget. */ 177 if (A != B) { 178 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 179 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 180 } 181 PetscFunctionReturn(PETSC_SUCCESS); 182 } 183 184 /*@C 185 DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation. 186 187 Logically Collective 188 189 Input Parameters: 190 + dm - `DM` to associate callback with 191 . func - local function evaluation 192 - ctx - context for function evaluation 193 194 Level: intermediate 195 196 Notes: 197 `func` should set the essential boundary data for the local portion of the solution, as 198 well its time derivative (if it is not `NULL`). 199 200 Vectors are initialized to zero before this function, so it is only needed for non 201 homogeneous data. 202 203 This function is somewhat optional: boundary data could potentially be inserted by a function 204 passed to `DMTSSetIFunctionLocal()`. The use case for this function is for discretizations 205 with constraints (see `DMGetDefaultConstraints()`): this function inserts boundary values 206 before constraint interpolation. 207 208 .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()` 209 @*/ 210 PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx) 211 { 212 DMTS tdm; 213 DMTS_Local *dmlocalts; 214 215 PetscFunctionBegin; 216 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 217 PetscCall(DMGetDMTSWrite(dm, &tdm)); 218 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 219 220 dmlocalts->boundarylocal = func; 221 dmlocalts->boundarylocalctx = ctx; 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 /*@C 226 DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector 227 containing the local vector information PLUS ghost point information. It should compute a result for all local 228 elements and `DM` will automatically accumulate the overlapping values. 229 230 Logically Collective 231 232 Input Parameter: 233 . dm - `DM` to associate callback with 234 235 Output Parameters: 236 + func - local function evaluation 237 - ctx - context for function evaluation 238 239 Level: beginner 240 241 .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()` 242 @*/ 243 PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx) 244 { 245 DMTS tdm; 246 DMTS_Local *dmlocalts; 247 248 PetscFunctionBegin; 249 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 250 PetscCall(DMGetDMTS(dm, &tdm)); 251 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 252 if (func) { 253 PetscAssertPointer(func, 2); 254 *func = dmlocalts->ifunctionlocal; 255 } 256 if (ctx) { 257 PetscAssertPointer(ctx, 3); 258 *ctx = dmlocalts->ifunctionlocalctx; 259 } 260 PetscFunctionReturn(PETSC_SUCCESS); 261 } 262 263 /*@C 264 DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector 265 containing the local vector information PLUS ghost point information. It should compute a result for all local 266 elements and `DM` will automatically accumulate the overlapping values. 267 268 Logically Collective 269 270 Input Parameters: 271 + dm - `DM` to associate callback with 272 . func - local function evaluation 273 - ctx - context for function evaluation 274 275 Level: beginner 276 277 .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()` 278 @*/ 279 PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx) 280 { 281 DMTS tdm; 282 DMTS_Local *dmlocalts; 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 286 PetscCall(DMGetDMTSWrite(dm, &tdm)); 287 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 288 289 dmlocalts->ifunctionlocal = func; 290 dmlocalts->ifunctionlocalctx = ctx; 291 292 PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts)); 293 if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 294 PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts)); 295 } 296 PetscFunctionReturn(PETSC_SUCCESS); 297 } 298 299 /*@C 300 DMTSGetIJacobianLocal - get a local Jacobian evaluation function 301 302 Logically Collective 303 304 Input Parameter: 305 . dm - `DM` to associate callback with 306 307 Output Parameters: 308 + func - local Jacobian evaluation 309 - ctx - optional context for local Jacobian evaluation 310 311 Level: beginner 312 313 .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()` 314 @*/ 315 PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx) 316 { 317 DMTS tdm; 318 DMTS_Local *dmlocalts; 319 320 PetscFunctionBegin; 321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 322 PetscCall(DMGetDMTS(dm, &tdm)); 323 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 324 if (func) { 325 PetscAssertPointer(func, 2); 326 *func = dmlocalts->ijacobianlocal; 327 } 328 if (ctx) { 329 PetscAssertPointer(ctx, 3); 330 *ctx = dmlocalts->ijacobianlocalctx; 331 } 332 PetscFunctionReturn(PETSC_SUCCESS); 333 } 334 335 /*@C 336 DMTSSetIJacobianLocal - set a local Jacobian evaluation function 337 338 Logically Collective 339 340 Input Parameters: 341 + dm - `DM` to associate callback with 342 . func - local Jacobian evaluation 343 - ctx - optional context for local Jacobian evaluation 344 345 Level: beginner 346 347 .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()` 348 @*/ 349 PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx) 350 { 351 DMTS tdm; 352 DMTS_Local *dmlocalts; 353 354 PetscFunctionBegin; 355 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 356 PetscCall(DMGetDMTSWrite(dm, &tdm)); 357 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 358 359 dmlocalts->ijacobianlocal = func; 360 dmlocalts->ijacobianlocalctx = ctx; 361 362 PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts)); 363 PetscFunctionReturn(PETSC_SUCCESS); 364 } 365 366 /*@C 367 DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector 368 containing the local vector information PLUS ghost point information. It should compute a result for all local 369 elements and `DM` will automatically accumulate the overlapping values. 370 371 Logically Collective 372 373 Input Parameter: 374 . dm - `DM` to associate callback with 375 376 Output Parameters: 377 + func - local function evaluation 378 - ctx - context for function evaluation 379 380 Level: beginner 381 382 .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()` 383 @*/ 384 PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx) 385 { 386 DMTS tdm; 387 DMTS_Local *dmlocalts; 388 389 PetscFunctionBegin; 390 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 391 PetscCall(DMGetDMTS(dm, &tdm)); 392 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 393 if (func) { 394 PetscAssertPointer(func, 2); 395 *func = dmlocalts->rhsfunctionlocal; 396 } 397 if (ctx) { 398 PetscAssertPointer(ctx, 3); 399 *ctx = dmlocalts->rhsfunctionlocalctx; 400 } 401 PetscFunctionReturn(PETSC_SUCCESS); 402 } 403 404 /*@C 405 DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector 406 containing the local vector information PLUS ghost point information. It should compute a result for all local 407 elements and `DM` will automatically accumulate the overlapping values. 408 409 Logically Collective 410 411 Input Parameters: 412 + dm - `DM` to associate callback with 413 . func - local function evaluation 414 - ctx - context for function evaluation 415 416 Level: beginner 417 418 .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()` 419 @*/ 420 PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx) 421 { 422 DMTS tdm; 423 DMTS_Local *dmlocalts; 424 425 PetscFunctionBegin; 426 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 427 PetscCall(DMGetDMTSWrite(dm, &tdm)); 428 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 429 430 dmlocalts->rhsfunctionlocal = func; 431 dmlocalts->rhsfunctionlocalctx = ctx; 432 433 PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts)); 434 PetscFunctionReturn(PETSC_SUCCESS); 435 } 436 437 /*@C 438 DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context. 439 440 Collective 441 442 Input Parameter: 443 . dm - `DM` providing the mass matrix 444 445 Level: developer 446 447 Note: 448 The idea here is that an explicit system can be given a mass matrix, based on the `DM`, which is inverted on the RHS at each step. 449 450 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS` 451 @*/ 452 PetscErrorCode DMTSCreateRHSMassMatrix(DM dm) 453 { 454 DMTS tdm; 455 DMTS_Local *dmlocalts; 456 const char *prefix; 457 458 PetscFunctionBegin; 459 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 460 PetscCall(DMGetDMTSWrite(dm, &tdm)); 461 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 462 PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass)); 463 PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass)); 464 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 465 PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix)); 466 PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_")); 467 PetscCall(KSPSetFromOptions(dmlocalts->kspmass)); 468 PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass)); 469 PetscFunctionReturn(PETSC_SUCCESS); 470 } 471 472 /*@C 473 DMTSCreateRHSMassMatrixLumped - This creates the lumped mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context. 474 475 Collective 476 477 Input Parameter: 478 . dm - `DM` providing the mass matrix 479 480 Level: developer 481 482 Note: 483 The idea here is that an explicit system can be given a mass matrix, based on the `DM`, which is inverted on the RHS at each step. 484 Since the matrix is lumped, inversion is trivial. 485 486 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS` 487 @*/ 488 PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm) 489 { 490 DMTS tdm; 491 DMTS_Local *dmlocalts; 492 493 PetscFunctionBegin; 494 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 495 PetscCall(DMGetDMTSWrite(dm, &tdm)); 496 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 497 PetscCall(DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv)); 498 PetscCall(VecReciprocal(dmlocalts->lumpedmassinv)); 499 PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view")); 500 PetscFunctionReturn(PETSC_SUCCESS); 501 } 502 503 /*@C 504 DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist. 505 506 Logically Collective 507 508 Input Parameter: 509 . dm - `DM` providing the mass matrix 510 511 Level: developer 512 513 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMTS` 514 @*/ 515 PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm) 516 { 517 DMTS tdm; 518 DMTS_Local *dmlocalts; 519 520 PetscFunctionBegin; 521 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 522 PetscCall(DMGetDMTSWrite(dm, &tdm)); 523 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 524 PetscCall(VecDestroy(&dmlocalts->lumpedmassinv)); 525 PetscCall(MatDestroy(&dmlocalts->mass)); 526 PetscCall(KSPDestroy(&dmlocalts->kspmass)); 527 PetscFunctionReturn(PETSC_SUCCESS); 528 } 529