Actual source code: dmlocalts.c

  1: #include <petsc/private/dmimpl.h>
  2: #include <petsc/private/tsimpl.h>

  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;

 18: static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
 19: {
 20:   PetscFree(tdm->data);
 21:   return 0;
 22: }

 24: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
 25: {
 26:   PetscNew((DMTS_Local **)&tdm->data);
 27:   if (oldtdm->data) PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local));
 28:   return 0;
 29: }

 31: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
 32: {
 33:   *dmlocalts = NULL;
 34:   if (!tdm->data) {
 35:     PetscNew((DMTS_Local **)&tdm->data);

 37:     tdm->ops->destroy   = DMTSDestroy_DMLocal;
 38:     tdm->ops->duplicate = DMTSDuplicate_DMLocal;
 39:   }
 40:   *dmlocalts = (DMTS_Local *)tdm->data;
 41:   return 0;
 42: }

 44: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx)
 45: {
 46:   DM          dm;
 47:   Vec         locX, locX_t, locF;
 48:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

 54:   TSGetDM(ts, &dm);
 55:   DMGetLocalVector(dm, &locX);
 56:   DMGetLocalVector(dm, &locX_t);
 57:   DMGetLocalVector(dm, &locF);
 58:   VecZeroEntries(locX);
 59:   VecZeroEntries(locX_t);
 60:   if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx);
 61:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
 62:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
 63:   DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
 64:   DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
 65:   VecZeroEntries(locF);
 66:   CHKMEMQ;
 67:   (*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx);
 68:   CHKMEMQ;
 69:   VecZeroEntries(F);
 70:   DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
 71:   DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
 72:   DMRestoreLocalVector(dm, &locX);
 73:   DMRestoreLocalVector(dm, &locX_t);
 74:   DMRestoreLocalVector(dm, &locF);
 75:   return 0;
 76: }

 78: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
 79: {
 80:   DM          dm;
 81:   Vec         locX, locF;
 82:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

 87:   TSGetDM(ts, &dm);
 88:   DMGetLocalVector(dm, &locX);
 89:   DMGetLocalVector(dm, &locF);
 90:   VecZeroEntries(locX);
 91:   if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx);
 92:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
 93:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
 94:   VecZeroEntries(locF);
 95:   CHKMEMQ;
 96:   (*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx);
 97:   CHKMEMQ;
 98:   VecZeroEntries(F);
 99:   DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
100:   DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
101:   if (dmlocalts->lumpedmassinv) {
102:     VecPointwiseMult(F, dmlocalts->lumpedmassinv, F);
103:   } else if (dmlocalts->kspmass) {
104:     Vec tmp;

106:     DMGetGlobalVector(dm, &tmp);
107:     KSPSolve(dmlocalts->kspmass, F, tmp);
108:     VecCopy(tmp, F);
109:     DMRestoreGlobalVector(dm, &tmp);
110:   }
111:   DMRestoreLocalVector(dm, &locX);
112:   DMRestoreLocalVector(dm, &locF);
113:   return 0;
114: }

116: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx)
117: {
118:   DM          dm;
119:   Vec         locX, locX_t;
120:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

122:   TSGetDM(ts, &dm);
123:   if (dmlocalts->ijacobianlocal) {
124:     DMGetLocalVector(dm, &locX);
125:     DMGetLocalVector(dm, &locX_t);
126:     VecZeroEntries(locX);
127:     VecZeroEntries(locX_t);
128:     if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx);
129:     DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
130:     DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
131:     DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
132:     DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
133:     CHKMEMQ;
134:     (*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx);
135:     CHKMEMQ;
136:     DMRestoreLocalVector(dm, &locX);
137:     DMRestoreLocalVector(dm, &locX_t);
138:   } else {
139:     MatFDColoring fdcoloring;
140:     PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring);
141:     if (!fdcoloring) {
142:       ISColoring coloring;

144:       DMCreateColoring(dm, dm->coloringtype, &coloring);
145:       MatFDColoringCreate(B, coloring, &fdcoloring);
146:       ISColoringDestroy(&coloring);
147:       switch (dm->coloringtype) {
148:       case IS_COLORING_GLOBAL:
149:         MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))TSComputeIFunction_DMLocal, dmlocalts);
150:         break;
151:       default:
152:         SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
153:       }
154:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix);
155:       MatFDColoringSetFromOptions(fdcoloring);
156:       MatFDColoringSetUp(B, coloring, fdcoloring);
157:       PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring);
158:       PetscObjectDereference((PetscObject)fdcoloring);

160:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
161:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
162:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
163:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
164:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
165:        */
166:       PetscObjectDereference((PetscObject)dm);
167:     }
168:     MatFDColoringApply(B, fdcoloring, X, ts);
169:   }
170:   /* This will be redundant if the user called both, but it's too common to forget. */
171:   if (A != B) {
172:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
173:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
174:   }
175:   return 0;
176: }

178: /*@C
179:   DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
180:     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).
181:     Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.

183:   Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
184:   `DMTSSetIFunctionLocal()`.  The use case for this function is for discretizations with constraints (see
185:   `DMGetDefaultConstraints()`): this function inserts boundary values before constraint interpolation.

187:   Logically Collective

189:   Input Parameters:
190: + dm   - `DM` to associate callback with
191: . func - local function evaluation
192: - ctx  - context for function evaluation

194:   Level: intermediate

196: .seealso: [](chapter_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
197: @*/
198: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
199: {
200:   DMTS        tdm;
201:   DMTS_Local *dmlocalts;

204:   DMGetDMTSWrite(dm, &tdm);
205:   DMLocalTSGetContext(dm, tdm, &dmlocalts);

207:   dmlocalts->boundarylocal    = func;
208:   dmlocalts->boundarylocalctx = ctx;

210:   return 0;
211: }

213: /*@C
214:   DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
215:       containing the local vector information PLUS ghost point information. It should compute a result for all local
216:       elements and `DM` will automatically accumulate the overlapping values.

218:   Logically Collective

220:   Input Parameter:
221: . dm   - `DM` to associate callback with

223:   Output Parameters:
224: + func - local function evaluation
225: - ctx  - context for function evaluation

227:   Level: beginner

229: .seealso: [](chapter_ts), `DM`, `DMTSSetIFunctionLocal(()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
230: @*/
231: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
232: {
233:   DMTS           tdm;
234:   DMTS_Local    *dmlocalts;

238:   DMGetDMTS(dm, &tdm);
239:   CHKERRQ(ierr);
240:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
241:   CHKERRQ(ierr);
242:   if (func) {
244:     *func = dmlocalts->ifunctionlocal;
245:   }
246:   if (ctx) {
248:     *ctx = dmlocalts->ifunctionlocalctx;
249:   }
250:   return 0;
251: }

253: /*@C
254:   DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
255:       containing the local vector information PLUS ghost point information. It should compute a result for all local
256:       elements and `DM` will automatically accumulate the overlapping values.

258:   Logically Collective

260:   Input Parameters:
261: + dm   - `DM` to associate callback with
262: . func - local function evaluation
263: - ctx  - context for function evaluation

265:   Level: beginner

267: .seealso: [](chapter_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
268: @*/
269: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
270: {
271:   DMTS        tdm;
272:   DMTS_Local *dmlocalts;

275:   DMGetDMTSWrite(dm, &tdm);
276:   DMLocalTSGetContext(dm, tdm, &dmlocalts);

278:   dmlocalts->ifunctionlocal    = func;
279:   dmlocalts->ifunctionlocalctx = ctx;

281:   DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts);
282:   if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
283:     DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
284:   }
285:   return 0;
286: }

288: /*@C
289:   DMTSGetIJacobianLocal - get a local Jacobian evaluation function

291:   Logically Collective

293:   Input Parameter:
294: . dm - `DM` to associate callback with

296:   Output Parameters:
297: + func - local Jacobian evaluation
298: - ctx - optional context for local Jacobian evaluation

300:   Level: beginner

302: .seealso: [](chapter_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
303: @*/
304: PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx)
305: {
306:   DMTS           tdm;
307:   DMTS_Local    *dmlocalts;

311:   DMGetDMTS(dm, &tdm);
312:   CHKERRQ(ierr);
313:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
314:   CHKERRQ(ierr);
315:   if (func) {
317:     *func = dmlocalts->ijacobianlocal;
318:   }
319:   if (ctx) {
321:     *ctx = dmlocalts->ijacobianlocalctx;
322:   }
323:   return 0;
324: }

326: /*@C
327:   DMTSSetIJacobianLocal - set a local Jacobian evaluation function

329:   Logically Collective

331:   Input Parameters:
332: + dm - `DM` to associate callback with
333: . func - local Jacobian evaluation
334: - ctx - optional context for local Jacobian evaluation

336:   Level: beginner

338: .seealso: [](chapter_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
339: @*/
340: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
341: {
342:   DMTS        tdm;
343:   DMTS_Local *dmlocalts;

346:   DMGetDMTSWrite(dm, &tdm);
347:   DMLocalTSGetContext(dm, tdm, &dmlocalts);

349:   dmlocalts->ijacobianlocal    = func;
350:   dmlocalts->ijacobianlocalctx = ctx;

352:   DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
353:   return 0;
354: }

356: /*@C
357:   DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
358:       containing the local vector information PLUS ghost point information. It should compute a result for all local
359:       elements and `DM` will automatically accumulate the overlapping values.

361:   Logically Collective

363:   Input Parameter:
364: . dm   - `DM` to associate callback with

366:   Output Parameters:
367: + func - local function evaluation
368: - ctx  - context for function evaluation

370:   Level: beginner

372: .seealso: [](chapter_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
373: @*/
374: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
375: {
376:   DMTS           tdm;
377:   DMTS_Local    *dmlocalts;

381:   DMGetDMTS(dm, &tdm);
382:   CHKERRQ(ierr);
383:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
384:   CHKERRQ(ierr);
385:   if (func) {
387:     *func = dmlocalts->rhsfunctionlocal;
388:   }
389:   if (ctx) {
391:     *ctx = dmlocalts->rhsfunctionlocalctx;
392:   }
393:   return 0;
394: }

396: /*@C
397:   DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
398:       containing the local vector information PLUS ghost point information. It should compute a result for all local
399:       elements and `DM` will automatically accumulate the overlapping values.

401:   Logically Collective

403:   Input Parameters:
404: + dm   - `DM` to associate callback with
405: . func - local function evaluation
406: - ctx  - context for function evaluation

408:   Level: beginner

410: .seealso: [](chapter_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
411: @*/
412: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
413: {
414:   DMTS        tdm;
415:   DMTS_Local *dmlocalts;

418:   DMGetDMTSWrite(dm, &tdm);
419:   DMLocalTSGetContext(dm, tdm, &dmlocalts);

421:   dmlocalts->rhsfunctionlocal    = func;
422:   dmlocalts->rhsfunctionlocalctx = ctx;

424:   DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts);
425:   return 0;
426: }

428: /*@C
429:   DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.

431:   Collective

433:   Input Parameters:
434: . dm   - `DM` providing the mass matrix

436:   Level: developer

438:   Note:
439:   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.

441: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
442: @*/
443: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
444: {
445:   DMTS        tdm;
446:   DMTS_Local *dmlocalts;
447:   const char *prefix;

450:   DMGetDMTSWrite(dm, &tdm);
451:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
452:   DMCreateMassMatrix(dm, dm, &dmlocalts->mass);
453:   KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass);
454:   PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
455:   KSPSetOptionsPrefix(dmlocalts->kspmass, prefix);
456:   KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_");
457:   KSPSetFromOptions(dmlocalts->kspmass);
458:   KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass);
459:   return 0;
460: }

462: /*@C
463:   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.

465:   Collective

467:   Input Parameters:
468: . dm   - `DM` providing the mass matrix

470:   Level: developer

472:   Note:
473:   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.
474:   Since the matrix is lumped, inversion is trivial.

476: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
477: @*/
478: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
479: {
480:   DMTS        tdm;
481:   DMTS_Local *dmlocalts;

484:   DMGetDMTSWrite(dm, &tdm);
485:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
486:   DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv);
487:   VecReciprocal(dmlocalts->lumpedmassinv);
488:   VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view");
489:   return 0;
490: }

492: /*@C
493:   DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.

495:   Logically Collective

497:   Input Parameters:
498: . dm   - `DM` providing the mass matrix

500:   Level: developer

502: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
503: @*/
504: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
505: {
506:   DMTS        tdm;
507:   DMTS_Local *dmlocalts;

510:   DMGetDMTSWrite(dm, &tdm);
511:   DMLocalTSGetContext(dm, tdm, &dmlocalts);
512:   VecDestroy(&dmlocalts->lumpedmassinv);
513:   MatDestroy(&dmlocalts->mass);
514:   KSPDestroy(&dmlocalts->kspmass);
515:   return 0;
516: }