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(0); 23 } 24 25 static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm) 26 { 27 PetscFunctionBegin; 28 PetscCall(PetscNewLog(tdm, (DMTS_Local **) &tdm->data)); 29 if (oldtdm->data) PetscCall(PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local))); 30 PetscFunctionReturn(0); 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(PetscNewLog(dm, (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(0); 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(0); 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(0); 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: SETERRQ(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]); 158 } 159 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) fdcoloring, ((PetscObject) dm)->prefix)); 160 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 161 PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring)); 162 PetscCall(PetscObjectCompose((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject) fdcoloring)); 163 PetscCall(PetscObjectDereference((PetscObject) fdcoloring)); 164 165 /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call 166 * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the 167 * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually 168 * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be 169 * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed. 170 */ 171 PetscCall(PetscObjectDereference((PetscObject) dm)); 172 } 173 PetscCall(MatFDColoringApply(B, fdcoloring, X, ts)); 174 } 175 /* This will be redundant if the user called both, but it's too common to forget. */ 176 if (A != B) { 177 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 178 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 179 } 180 PetscFunctionReturn(0); 181 } 182 183 /*@C 184 DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation. 185 It should set the essential boundary data for the local portion of the solution X, as well its time derivative X_t (if it is not NULL). 186 Vectors are initialized to zero before this function, so it is only needed for non homogeneous data. 187 188 Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to 189 DMTSSetIFunctionLocal(). The use case for this function is for discretizations with constraints (see 190 DMGetDefaultConstraints()): this function inserts boundary values before constraint interpolation. 191 192 Logically Collective 193 194 Input Parameters: 195 + dm - DM to associate callback with 196 . func - local function evaluation 197 - ctx - context for function evaluation 198 199 Level: intermediate 200 201 .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal() 202 @*/ 203 PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx) 204 { 205 DMTS tdm; 206 DMTS_Local *dmlocalts; 207 208 PetscFunctionBegin; 209 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 210 PetscCall(DMGetDMTSWrite(dm, &tdm)); 211 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 212 213 dmlocalts->boundarylocal = func; 214 dmlocalts->boundarylocalctx = ctx; 215 216 PetscFunctionReturn(0); 217 } 218 219 /*@C 220 DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector 221 containing the local vector information PLUS ghost point information. It should compute a result for all local 222 elements and DMTS will automatically accumulate the overlapping values. 223 224 Logically Collective 225 226 Input Parameters: 227 + dm - DM to associate callback with 228 . func - local function evaluation 229 - ctx - context for function evaluation 230 231 Level: beginner 232 233 .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal() 234 @*/ 235 PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx) 236 { 237 DMTS tdm; 238 DMTS_Local *dmlocalts; 239 240 PetscFunctionBegin; 241 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 242 PetscCall(DMGetDMTSWrite(dm, &tdm)); 243 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 244 245 dmlocalts->ifunctionlocal = func; 246 dmlocalts->ifunctionlocalctx = ctx; 247 248 PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts)); 249 if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */ 250 PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts)); 251 } 252 PetscFunctionReturn(0); 253 } 254 255 /*@C 256 DMTSSetIJacobianLocal - set a local Jacobian evaluation function 257 258 Logically Collective 259 260 Input Parameters: 261 + dm - DM to associate callback with 262 . func - local Jacobian evaluation 263 - ctx - optional context for local Jacobian evaluation 264 265 Level: beginner 266 267 .seealso: DMTSSetIFunctionLocal(), DMTSSetIJacobian(), DMTSSetIFunction() 268 @*/ 269 PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx) 270 { 271 DMTS tdm; 272 DMTS_Local *dmlocalts; 273 274 PetscFunctionBegin; 275 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 276 PetscCall(DMGetDMTSWrite(dm, &tdm)); 277 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 278 279 dmlocalts->ijacobianlocal = func; 280 dmlocalts->ijacobianlocalctx = ctx; 281 282 PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts)); 283 PetscFunctionReturn(0); 284 } 285 286 /*@C 287 DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector 288 containing the local vector information PLUS ghost point information. It should compute a result for all local 289 elements and DMTS will automatically accumulate the overlapping values. 290 291 Logically Collective 292 293 Input Parameters: 294 + dm - DM to associate callback with 295 . func - local function evaluation 296 - ctx - context for function evaluation 297 298 Level: beginner 299 300 .seealso: DMTSSetRHSFunction(), DMTSSetIFunction(), DMTSSetIJacobianLocal() 301 @*/ 302 PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx) 303 { 304 DMTS tdm; 305 DMTS_Local *dmlocalts; 306 307 PetscFunctionBegin; 308 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 309 PetscCall(DMGetDMTSWrite(dm, &tdm)); 310 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 311 312 dmlocalts->rhsfunctionlocal = func; 313 dmlocalts->rhsfunctionlocalctx = ctx; 314 315 PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts)); 316 PetscFunctionReturn(0); 317 } 318 319 /*@C 320 DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context. 321 322 Collective on dm 323 324 Input Parameters: 325 . dm - DM providing the mass matrix 326 327 Note: 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. 328 329 Level: developer 330 331 .seealso: DMTSCreateRHSMassMatrixLumped(), DMTSDestroyRHSMassMatrix(), DMCreateMassMatrix(), DMTS 332 @*/ 333 PetscErrorCode DMTSCreateRHSMassMatrix(DM dm) 334 { 335 DMTS tdm; 336 DMTS_Local *dmlocalts; 337 const char *prefix; 338 339 PetscFunctionBegin; 340 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 341 PetscCall(DMGetDMTSWrite(dm, &tdm)); 342 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 343 PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass)); 344 PetscCall(KSPCreate(PetscObjectComm((PetscObject) dm), &dmlocalts->kspmass)); 345 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 346 PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix)); 347 PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_")); 348 PetscCall(KSPSetFromOptions(dmlocalts->kspmass)); 349 PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass)); 350 PetscFunctionReturn(0); 351 } 352 353 /*@C 354 DMTSCreateRHSMassMatrixLumped - This creates the lumped mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context. 355 356 Collective on dm 357 358 Input Parameters: 359 . dm - DM providing the mass matrix 360 361 Note: 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. 362 Since the matrix is lumped, inversion is trivial. 363 364 Level: developer 365 366 .seealso: DMTSCreateRHSMassMatrix(), DMTSDestroyRHSMassMatrix(), DMCreateMassMatrix(), DMTS 367 @*/ 368 PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm) 369 { 370 DMTS tdm; 371 DMTS_Local *dmlocalts; 372 373 PetscFunctionBegin; 374 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 375 PetscCall(DMGetDMTSWrite(dm, &tdm)); 376 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 377 PetscCall(DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv)); 378 PetscCall(VecReciprocal(dmlocalts->lumpedmassinv)); 379 PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view")); 380 PetscFunctionReturn(0); 381 } 382 383 /*@C 384 DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the DMTS context, if they exist. 385 386 Logically Collective 387 388 Input Parameters: 389 . dm - DM providing the mass matrix 390 391 Level: developer 392 393 .seealso: DMTSCreateRHSMassMatrixLumped(), DMCreateMassMatrix(), DMCreateMassMatrix(), DMTS 394 @*/ 395 PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm) 396 { 397 DMTS tdm; 398 DMTS_Local *dmlocalts; 399 400 PetscFunctionBegin; 401 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 402 PetscCall(DMGetDMTSWrite(dm, &tdm)); 403 PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts)); 404 PetscCall(VecDestroy(&dmlocalts->lumpedmassinv)); 405 PetscCall(MatDestroy(&dmlocalts->mass)); 406 PetscCall(KSPDestroy(&dmlocalts->kspmass)); 407 PetscFunctionReturn(0); 408 } 409