Actual source code: dmadapt.c

  1: #include <petscdmadaptor.h>
  2: #include <petscdmplex.h>
  3: #include <petscdmforest.h>
  4: #include <petscds.h>
  5: #include <petscblaslapack.h>

  7: #include <petsc/private/dmadaptorimpl.h>
  8: #include <petsc/private/dmpleximpl.h>
  9: #include <petsc/private/petscfeimpl.h>

 11: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor, PetscInt, PetscInt, const PetscScalar *, const PetscScalar *, const PetscFVCellGeom *, PetscReal *, void *);

 13: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
 14: {
 15:   PetscFunctionBeginUser;
 16:   PetscCall(DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: /*@
 21:   DMAdaptorCreate - Create a `DMAdaptor` object. Its purpose is to construct a adaptation `DMLabel` or metric Vec that can be used to modify the `DM`.

 23:   Collective

 25:   Input Parameter:
 26: . comm - The communicator for the `DMAdaptor` object

 28:   Output Parameter:
 29: . adaptor   - The `DMAdaptor` object

 31:   Level: beginner

 33: .seealso: `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`
 34: @*/
 35: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
 36: {
 37:   VecTaggerBox refineBox, coarsenBox;

 39:   PetscFunctionBegin;
 41:   PetscCall(PetscSysInitializePackage());
 42:   PetscCall(PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView));

 44:   (*adaptor)->monitor                    = PETSC_FALSE;
 45:   (*adaptor)->adaptCriterion             = DM_ADAPTATION_NONE;
 46:   (*adaptor)->numSeq                     = 1;
 47:   (*adaptor)->Nadapt                     = -1;
 48:   (*adaptor)->refinementFactor           = 2.0;
 49:   (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
 50:   refineBox.min = refineBox.max = PETSC_MAX_REAL;
 51:   PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->refineTag));
 52:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->refineTag, "refine_"));
 53:   PetscCall(VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE));
 54:   PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox));
 55:   coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
 56:   PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->coarsenTag));
 57:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->coarsenTag, "coarsen_"));
 58:   PetscCall(VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE));
 59:   PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox));
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: /*@
 64:   DMAdaptorDestroy - Destroys a `DMAdaptor` object

 66:   Collective

 68:   Input Parameter:
 69: . adaptor - The `DMAdaptor` object

 71:   Level: beginner

 73: .seealso: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
 74: @*/
 75: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
 76: {
 77:   PetscFunctionBegin;
 78:   if (!*adaptor) PetscFunctionReturn(PETSC_SUCCESS);
 80:   if (--((PetscObject)(*adaptor))->refct > 0) {
 81:     *adaptor = NULL;
 82:     PetscFunctionReturn(PETSC_SUCCESS);
 83:   }
 84:   PetscCall(VecTaggerDestroy(&(*adaptor)->refineTag));
 85:   PetscCall(VecTaggerDestroy(&(*adaptor)->coarsenTag));
 86:   PetscCall(PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx));
 87:   PetscCall(PetscHeaderDestroy(adaptor));
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: /*@
 92:   DMAdaptorSetFromOptions - Sets properties of a `DMAdaptor` object from the options database

 94:   Collective

 96:   Input Parameter:
 97: . adaptor - The `DMAdaptor` object

 99:   Options Database Keys:
100: + -adaptor_monitor <bool>        - Monitor the adaptation process
101: . -adaptor_sequence_num <num>    - Number of adaptations to generate an optimal grid
102: . -adaptor_target_num <num>      - Set the target number of vertices N_adapt, -1 for automatic determination
103: - -adaptor_refinement_factor <r> - Set r such that N_adapt = r^dim N_orig

105:   Level: beginner

107: .seealso: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
108: @*/
109: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
110: {
111:   PetscFunctionBegin;
112:   PetscOptionsBegin(PetscObjectComm((PetscObject)adaptor), "", "DM Adaptor Options", "DMAdaptor");
113:   PetscCall(PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL));
114:   PetscCall(PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL));
115:   PetscCall(PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL));
116:   PetscCall(PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL));
117:   PetscOptionsEnd();
118:   PetscCall(VecTaggerSetFromOptions(adaptor->refineTag));
119:   PetscCall(VecTaggerSetFromOptions(adaptor->coarsenTag));
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: /*@
124:    DMAdaptorView - Views a `DMAdaptor` object

126:    Collective

128:    Input Parameters:
129: +  adaptor - The `DMAdaptor` object
130: -  viewer - The `PetscViewer` object

132:   Level: beginner

134: .seealso: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
135: @*/
136: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
137: {
138:   PetscFunctionBegin;
139:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer));
140:   PetscCall(PetscViewerASCIIPrintf(viewer, "DM Adaptor\n"));
141:   PetscCall(PetscViewerASCIIPrintf(viewer, "  sequence length: %" PetscInt_FMT "\n", adaptor->numSeq));
142:   PetscCall(VecTaggerView(adaptor->refineTag, viewer));
143:   PetscCall(VecTaggerView(adaptor->coarsenTag, viewer));
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

147: /*@
148:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

150:   Not Collective

152:   Input Parameter:
153: . adaptor   - The `DMAdaptor` object

155:   Output Parameter:
156: . snes - The solver

158:   Level: intermediate

160: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
161: @*/
162: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
163: {
164:   PetscFunctionBegin;
167:   *snes = adaptor->snes;
168:   PetscFunctionReturn(PETSC_SUCCESS);
169: }

171: /*@
172:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

174:   Not Collective

176:   Input Parameters:
177: + adaptor   - The `DMAdaptor` object
178: - snes - The solver

180:   Level: intermediate

182:   Note:
183:   The solver MUST have an attached `DM`/`PetscDS`, so that we know the exact solution

185: .seealso: `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
186: @*/
187: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
188: {
189:   PetscFunctionBegin;
192:   adaptor->snes = snes;
193:   PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
194:   PetscFunctionReturn(PETSC_SUCCESS);
195: }

197: /*@
198:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter

200:   Not Collective

202:   Input Parameter:
203: . adaptor - The `DMAdaptor` object

205:   Output Parameter:
206: . num - The number of adaptations

208:   Level: intermediate

210: .seealso: `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
211: @*/
212: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
213: {
214:   PetscFunctionBegin;
217:   *num = adaptor->numSeq;
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /*@
222:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

224:   Not Collective

226:   Input Parameters:
227: + adaptor - The `DMAdaptor` object
228: - num - The number of adaptations

230:   Level: intermediate

232: .seealso: `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
233: @*/
234: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
235: {
236:   PetscFunctionBegin;
238:   adaptor->numSeq = num;
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: /*@
243:   DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity

245:   Collective

247:   Input Parameter:
248: . adaptor - The `DMAdaptor` object

250:   Level: beginner

252: .seealso: `DMAdaptorCreate()`, `DMAdaptorAdapt()`
253: @*/
254: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
255: {
256:   PetscDS  prob;
257:   PetscInt Nf, f;

259:   PetscFunctionBegin;
260:   PetscCall(DMGetDS(adaptor->idm, &prob));
261:   PetscCall(VecTaggerSetUp(adaptor->refineTag));
262:   PetscCall(VecTaggerSetUp(adaptor->coarsenTag));
263:   PetscCall(PetscDSGetNumFields(prob, &Nf));
264:   PetscCall(PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx));
265:   for (f = 0; f < Nf; ++f) {
266:     PetscCall(PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]));
267:     /* TODO Have a flag that forces projection rather than using the exact solution */
268:     if (adaptor->exactSol[0]) PetscCall(DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private));
269:   }
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
274: {
275:   PetscFunctionBegin;
276:   *tfunc = adaptor->ops->transfersolution;
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
281: {
282:   PetscFunctionBegin;
283:   adaptor->ops->transfersolution = tfunc;
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
288: {
289:   DM           plex;
290:   PetscDS      prob;
291:   PetscObject  obj;
292:   PetscClassId id;
293:   PetscBool    isForest;

295:   PetscFunctionBegin;
296:   PetscCall(DMConvert(adaptor->idm, DMPLEX, &plex));
297:   PetscCall(DMGetDS(adaptor->idm, &prob));
298:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
299:   PetscCall(PetscObjectGetClassId(obj, &id));
300:   PetscCall(DMIsForest(adaptor->idm, &isForest));
301:   if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
302:     if (isForest) {
303:       adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
304:     }
305: #if defined(PETSC_HAVE_PRAGMATIC)
306:     else {
307:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
308:     }
309: #elif defined(PETSC_HAVE_MMG)
310:     else {
311:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
312:     }
313: #elif defined(PETSC_HAVE_PARMMG)
314:     else {
315:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
316:     }
317: #else
318:     else {
319:       adaptor->adaptCriterion = DM_ADAPTATION_REFINE;
320:     }
321: #endif
322:   }
323:   if (id == PETSCFV_CLASSID) {
324:     adaptor->femType = PETSC_FALSE;
325:   } else {
326:     adaptor->femType = PETSC_TRUE;
327:   }
328:   if (adaptor->femType) {
329:     /* Compute local solution bc */
330:     PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
331:   } else {
332:     PetscFV      fvm = (PetscFV)obj;
333:     PetscLimiter noneLimiter;
334:     Vec          grad;

336:     PetscCall(PetscFVGetComputeGradients(fvm, &adaptor->computeGradient));
337:     PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
338:     /* Use no limiting when reconstructing gradients for adaptivity */
339:     PetscCall(PetscFVGetLimiter(fvm, &adaptor->limiter));
340:     PetscCall(PetscObjectReference((PetscObject)adaptor->limiter));
341:     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
342:     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
343:     PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
344:     /* Get FVM data */
345:     PetscCall(DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM));
346:     PetscCall(VecGetDM(adaptor->cellGeom, &adaptor->cellDM));
347:     PetscCall(VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
348:     /* Compute local solution bc */
349:     PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
350:     /* Compute gradients */
351:     PetscCall(DMCreateGlobalVector(adaptor->gradDM, &grad));
352:     PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
353:     PetscCall(DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad));
354:     PetscCall(DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
355:     PetscCall(DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
356:     PetscCall(VecDestroy(&grad));
357:     PetscCall(VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
358:   }
359:   PetscCall(DMDestroy(&plex));
360:   PetscFunctionReturn(PETSC_SUCCESS);
361: }

363: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
364: {
365:   PetscReal time = 0.0;
366:   Mat       interp;
367:   void     *ctx;

369:   PetscFunctionBegin;
370:   PetscCall(DMGetApplicationContext(dm, &ctx));
371:   if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
372:   else {
373:     switch (adaptor->adaptCriterion) {
374:     case DM_ADAPTATION_LABEL:
375:       PetscCall(DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time));
376:       break;
377:     case DM_ADAPTATION_REFINE:
378:     case DM_ADAPTATION_METRIC:
379:       PetscCall(DMCreateInterpolation(dm, adm, &interp, NULL));
380:       PetscCall(MatInterpolate(interp, x, ax));
381:       PetscCall(DMInterpolate(dm, interp, adm));
382:       PetscCall(MatDestroy(&interp));
383:       break;
384:     default:
385:       SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
386:     }
387:   }
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
392: {
393:   PetscDS      prob;
394:   PetscObject  obj;
395:   PetscClassId id;

397:   PetscFunctionBegin;
398:   PetscCall(DMGetDS(adaptor->idm, &prob));
399:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
400:   PetscCall(PetscObjectGetClassId(obj, &id));
401:   if (id == PETSCFV_CLASSID) {
402:     PetscFV fvm = (PetscFV)obj;

404:     PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
405:     /* Restore original limiter */
406:     PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));

408:     PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
409:     PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
410:     PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
411:   }
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: /*
416:   DMAdaptorSimpleErrorIndicator - Use the integrated gradient as an error indicator

418:   Input Parameters:
419: + adaptor  - The `DMAdaptor` object
420: . dim      - The topological dimension
421: . cell     - The cell
422: . field    - The field integrated over the cell
423: . gradient - The gradient integrated over the cell
424: . cg       - A `PetscFVCellGeom` struct
425: - ctx      - A user context

427:   Output Parameter:
428: . errInd   - The error indicator

430:   Developer Note:
431:   Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API

433: .seealso: `DMAdaptor`, `DMAdaptorComputeErrorIndicator()`
434: */
435: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
436: {
437:   PetscReal err = 0.;
438:   PetscInt  c, d;

440:   PetscFunctionBeginHot;
441:   for (c = 0; c < Nc; c++) {
442:     for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
443:   }
444:   *errInd = cg->volume * err;
445:   PetscFunctionReturn(PETSC_SUCCESS);
446: }

448: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
449: {
450:   PetscDS         prob;
451:   PetscObject     obj;
452:   PetscClassId    id;
453:   void           *ctx;
454:   PetscQuadrature quad;
455:   PetscInt        dim, d, cdim, Nc;

457:   PetscFunctionBegin;
458:   *errInd = 0.;
459:   PetscCall(DMGetDimension(plex, &dim));
460:   PetscCall(DMGetCoordinateDim(plex, &cdim));
461:   PetscCall(DMGetApplicationContext(plex, &ctx));
462:   PetscCall(DMGetDS(plex, &prob));
463:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
464:   PetscCall(PetscObjectGetClassId(obj, &id));
465:   if (id == PETSCFV_CLASSID) {
466:     const PetscScalar *pointSols;
467:     const PetscScalar *pointSol;
468:     const PetscScalar *pointGrad;
469:     PetscFVCellGeom   *cg;

471:     PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
472:     PetscCall(VecGetArrayRead(locX, &pointSols));
473:     PetscCall(DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol));
474:     PetscCall(DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad));
475:     PetscCall(DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg));
476:     PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
477:     PetscCall(VecRestoreArrayRead(locX, &pointSols));
478:   } else {
479:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
480:     PetscFVCellGeom  cg;
481:     PetscFEGeom      fegeom;
482:     const PetscReal *quadWeights;
483:     PetscReal       *coords;
484:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

486:     fegeom.dim      = dim;
487:     fegeom.dimEmbed = cdim;
488:     PetscCall(PetscDSGetNumFields(prob, &Nf));
489:     PetscCall(PetscFEGetQuadrature((PetscFE)obj, &quad));
490:     PetscCall(DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x));
491:     PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
492:     PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
493:     PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights));
494:     PetscCall(PetscCalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ));
495:     PetscCall(PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad));
496:     PetscCall(DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
497:     PetscCall(DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL));
498:     PetscCall(PetscArrayzero(gradient, cdim * Nc));
499:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
500:       PetscInt qc = 0, q;

502:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
503:       PetscCall(PetscArrayzero(interpolant, Nc));
504:       PetscCall(PetscArrayzero(interpolantGrad, cdim * Nc));
505:       for (q = 0; q < Nq; ++q) {
506:         PetscCall(PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad));
507:         for (fc = 0; fc < Nc; ++fc) {
508:           const PetscReal wt = quadWeights[q * qNc + qc + fc];

510:           field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
511:           for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
512:         }
513:       }
514:       fieldOffset += Nb;
515:       qc += Nc;
516:     }
517:     PetscCall(PetscFree2(interpolant, interpolantGrad));
518:     PetscCall(DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x));
519:     for (fc = 0; fc < Nc; ++fc) {
520:       field[fc] /= cg.volume;
521:       for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
522:     }
523:     PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, field, gradient, &cg, errInd, ctx);
524:     PetscCall(PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
525:   }
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }

529: static void identityFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
530: {
531:   PetscInt i, j;

533:   for (i = 0; i < dim; ++i) {
534:     for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
535:   }
536: }

538: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
539: {
540:   PetscDS  prob;
541:   void    *ctx;
542:   MPI_Comm comm;
543:   PetscInt numAdapt = adaptor->numSeq, adaptIter;
544:   PetscInt dim, coordDim, numFields, cStart, cEnd, c;

546:   PetscFunctionBegin;
547:   PetscCall(DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view"));
548:   PetscCall(VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view"));
549:   PetscCall(PetscObjectGetComm((PetscObject)adaptor, &comm));
550:   PetscCall(DMGetDimension(adaptor->idm, &dim));
551:   PetscCall(DMGetCoordinateDim(adaptor->idm, &coordDim));
552:   PetscCall(DMGetApplicationContext(adaptor->idm, &ctx));
553:   PetscCall(DMGetDS(adaptor->idm, &prob));
554:   PetscCall(PetscDSGetNumFields(prob, &numFields));
555:   PetscCheck(numFields != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fields is zero!");

557:   /* Adapt until nothing changes */
558:   /* Adapt for a specified number of iterates */
559:   for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm)));
560:   for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
561:     PetscBool adapted = PETSC_FALSE;
562:     DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
563:     Vec       x       = adaptIter ? *ax : inx, locX, ox;

565:     PetscCall(DMGetLocalVector(dm, &locX));
566:     PetscCall(DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
567:     PetscCall(DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
568:     PetscCall(DMAdaptorPreAdapt(adaptor, locX));
569:     if (doSolve) {
570:       SNES snes;

572:       PetscCall(DMAdaptorGetSolver(adaptor, &snes));
573:       PetscCall(SNESSolve(snes, NULL, adaptIter ? *ax : x));
574:     }
575:     /* PetscCall(DMAdaptorMonitor(adaptor));
576:        Print iterate, memory used, DM, solution */
577:     switch (adaptor->adaptCriterion) {
578:     case DM_ADAPTATION_REFINE:
579:       PetscCall(DMRefine(dm, comm, &odm));
580:       PetscCheck(odm, comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
581:       adapted = PETSC_TRUE;
582:       break;
583:     case DM_ADAPTATION_LABEL: {
584:       /* Adapt DM
585:            Create local solution
586:            Reconstruct gradients (FVM) or solve adjoint equation (FEM)
587:            Produce cellwise error indicator */
588:       DM                 plex;
589:       DMLabel            adaptLabel;
590:       IS                 refineIS, coarsenIS;
591:       Vec                errVec;
592:       PetscScalar       *errArray;
593:       const PetscScalar *pointSols;
594:       PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
595:       PetscInt           nRefine, nCoarsen;

597:       PetscCall(DMConvert(dm, DMPLEX, &plex));
598:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
599:       PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));

601:       PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)adaptor), cEnd - cStart, PETSC_DETERMINE, &errVec));
602:       PetscCall(VecSetUp(errVec));
603:       PetscCall(VecGetArray(errVec, &errArray));
604:       PetscCall(VecGetArrayRead(locX, &pointSols));
605:       for (c = cStart; c < cEnd; ++c) {
606:         PetscReal errInd;

608:         PetscCall(DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd));
609:         errArray[c - cStart] = errInd;
610:         minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
611:         minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
612:       }
613:       PetscCall(VecRestoreArrayRead(locX, &pointSols));
614:       PetscCall(VecRestoreArray(errVec, &errArray));
615:       PetscCall(PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxIndGlobal));
616:       PetscCall(PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxIndGlobal[0], (double)minMaxIndGlobal[1]));
617:       /*     Compute IS from VecTagger */
618:       PetscCall(VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS, NULL));
619:       PetscCall(VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS, NULL));
620:       PetscCall(ISGetSize(refineIS, &nRefine));
621:       PetscCall(ISGetSize(coarsenIS, &nCoarsen));
622:       PetscCall(PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen));
623:       if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
624:       if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
625:       PetscCall(ISDestroy(&coarsenIS));
626:       PetscCall(ISDestroy(&refineIS));
627:       PetscCall(VecDestroy(&errVec));
628:       /*     Adapt DM from label */
629:       if (nRefine || nCoarsen) {
630:         PetscCall(DMAdaptLabel(dm, adaptLabel, &odm));
631:         adapted = PETSC_TRUE;
632:       }
633:       PetscCall(DMLabelDestroy(&adaptLabel));
634:       PetscCall(DMDestroy(&plex));
635:     } break;
636:     case DM_ADAPTATION_METRIC: {
637:       DM        dmGrad, dmHess, dmMetric, dmDet;
638:       Vec       xGrad, xHess, metric, determinant;
639:       PetscReal N;
640:       DMLabel   bdLabel = NULL, rgLabel = NULL;
641:       PetscBool higherOrder = PETSC_FALSE;
642:       PetscInt  Nd          = coordDim * coordDim, f, vStart, vEnd;
643:       void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

645:       PetscCall(PetscMalloc(1, &funcs));
646:       funcs[0] = identityFunc;

648:       /*     Setup finite element spaces */
649:       PetscCall(DMClone(dm, &dmGrad));
650:       PetscCall(DMClone(dm, &dmHess));
651:       PetscCheck(numFields <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple fields not yet considered"); // TODO
652:       for (f = 0; f < numFields; ++f) {
653:         PetscFE         fe, feGrad, feHess;
654:         PetscDualSpace  Q;
655:         PetscSpace      space;
656:         DM              K;
657:         PetscQuadrature q;
658:         PetscInt        Nc, qorder, p;
659:         const char     *prefix;

661:         PetscCall(PetscDSGetDiscretization(prob, f, (PetscObject *)&fe));
662:         PetscCall(PetscFEGetNumComponents(fe, &Nc));
663:         PetscCheck(Nc <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple components not yet considered"); // TODO
664:         PetscCall(PetscFEGetBasisSpace(fe, &space));
665:         PetscCall(PetscSpaceGetDegree(space, NULL, &p));
666:         if (p > 1) higherOrder = PETSC_TRUE;
667:         PetscCall(PetscFEGetDualSpace(fe, &Q));
668:         PetscCall(PetscDualSpaceGetDM(Q, &K));
669:         PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
670:         PetscCall(PetscFEGetQuadrature(fe, &q));
671:         PetscCall(PetscQuadratureGetOrder(q, &qorder));
672:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix));
673:         PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad));
674:         PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess));
675:         PetscCall(DMSetField(dmGrad, f, NULL, (PetscObject)feGrad));
676:         PetscCall(DMSetField(dmHess, f, NULL, (PetscObject)feHess));
677:         PetscCall(DMCreateDS(dmGrad));
678:         PetscCall(DMCreateDS(dmHess));
679:         PetscCall(PetscFEDestroy(&feGrad));
680:         PetscCall(PetscFEDestroy(&feHess));
681:       }
682:       /*     Compute vertexwise gradients from cellwise gradients */
683:       PetscCall(DMCreateLocalVector(dmGrad, &xGrad));
684:       PetscCall(VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view"));
685:       PetscCall(DMPlexComputeGradientClementInterpolant(dm, locX, xGrad));
686:       PetscCall(VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view"));
687:       /*     Compute vertexwise Hessians from cellwise Hessians */
688:       PetscCall(DMCreateLocalVector(dmHess, &xHess));
689:       PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess));
690:       PetscCall(VecViewFromOptions(xHess, NULL, "-adapt_hessian_view"));
691:       PetscCall(VecDestroy(&xGrad));
692:       PetscCall(DMDestroy(&dmGrad));
693:       /*     Compute L-p normalized metric */
694:       PetscCall(DMClone(dm, &dmMetric));
695:       N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
696:       if (adaptor->monitor) {
697:         PetscMPIInt rank, size;
698:         PetscCallMPI(MPI_Comm_rank(comm, &size));
699:         PetscCallMPI(MPI_Comm_rank(comm, &rank));
700:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] N_orig: %" PetscInt_FMT " N_adapt: %g\n", rank, vEnd - vStart, (double)N));
701:       }
702:       PetscCall(DMPlexMetricSetTargetComplexity(dmMetric, (PetscReal)N));
703:       if (higherOrder) {
704:         /*   Project Hessian into P1 space, if required */
705:         PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
706:         PetscCall(DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric));
707:         PetscCall(VecDestroy(&xHess));
708:         xHess = metric;
709:       }
710:       PetscCall(PetscFree(funcs));
711:       PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
712:       PetscCall(DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet));
713:       PetscCall(DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant));
714:       PetscCall(VecDestroy(&determinant));
715:       PetscCall(DMDestroy(&dmDet));
716:       PetscCall(VecDestroy(&xHess));
717:       PetscCall(DMDestroy(&dmHess));
718:       /*     Adapt DM from metric */
719:       PetscCall(DMGetLabel(dm, "marker", &bdLabel));
720:       PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm));
721:       adapted = PETSC_TRUE;
722:       /*     Cleanup */
723:       PetscCall(VecDestroy(&metric));
724:       PetscCall(DMDestroy(&dmMetric));
725:     } break;
726:     default:
727:       SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
728:     }
729:     PetscCall(DMAdaptorPostAdapt(adaptor));
730:     PetscCall(DMRestoreLocalVector(dm, &locX));
731:     /* If DM was adapted, replace objects and recreate solution */
732:     if (adapted) {
733:       const char *name;

735:       PetscCall(PetscObjectGetName((PetscObject)dm, &name));
736:       PetscCall(PetscObjectSetName((PetscObject)odm, name));
737:       /* Reconfigure solver */
738:       PetscCall(SNESReset(adaptor->snes));
739:       PetscCall(SNESSetDM(adaptor->snes, odm));
740:       PetscCall(DMAdaptorSetSolver(adaptor, adaptor->snes));
741:       PetscCall(DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx));
742:       PetscCall(SNESSetFromOptions(adaptor->snes));
743:       /* Transfer system */
744:       PetscCall(DMCopyDisc(dm, odm));
745:       /* Transfer solution */
746:       PetscCall(DMCreateGlobalVector(odm, &ox));
747:       PetscCall(PetscObjectGetName((PetscObject)x, &name));
748:       PetscCall(PetscObjectSetName((PetscObject)ox, name));
749:       PetscCall(DMAdaptorTransferSolution(adaptor, dm, x, odm, ox));
750:       /* Cleanup adaptivity info */
751:       if (adaptIter > 0) PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm)));
752:       PetscCall(DMForestSetAdaptivityForest(dm, NULL)); /* clear internal references to the previous dm */
753:       PetscCall(DMDestroy(&dm));
754:       PetscCall(VecDestroy(&x));
755:       *adm = odm;
756:       *ax  = ox;
757:     } else {
758:       *adm      = dm;
759:       *ax       = x;
760:       adaptIter = numAdapt;
761:     }
762:     if (adaptIter < numAdapt - 1) {
763:       PetscCall(DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view"));
764:       PetscCall(VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view"));
765:     }
766:   }
767:   PetscCall(DMViewFromOptions(*adm, NULL, "-dm_adapt_view"));
768:   PetscCall(VecViewFromOptions(*ax, NULL, "-sol_adapt_view"));
769:   PetscFunctionReturn(PETSC_SUCCESS);
770: }

772: /*@
773:   DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem

775:   Not Collective

777:   Input Parameters:
778: + adaptor  - The `DMAdaptor` object
779: . x        - The global approximate solution
780: - strategy - The adaptation strategy

782:   Output Parameters:
783: + adm - The adapted `DM`
784: - ax  - The adapted solution

786:   Options database Keys:
787: + -snes_adapt <strategy> - initial, sequential, multigrid
788: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
789: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
790: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement

792:   Level: intermediate

794:   Note:
795:   The available adaptation strategies are:
796: + * - Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
797: . * - Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
798: - * - Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid

800: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
801: @*/
802: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
803: {
804:   PetscFunctionBegin;
805:   switch (strategy) {
806:   case DM_ADAPTATION_INITIAL:
807:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
808:     break;
809:   case DM_ADAPTATION_SEQUENTIAL:
810:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
811:     break;
812:   default:
813:     SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
814:   }
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }