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: `DM`, `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`, `PetscConvEst`, `PetscConvEstCreate()`
 34: @*/
 35: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
 36: {
 37:   VecTaggerBox refineBox, coarsenBox;

 39:   PetscFunctionBegin;
 40:   PetscAssertPointer(adaptor, 2);
 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: `DM`, `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 values in 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: `DM`, `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: `DM`, `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: `DM`, `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
161: @*/
162: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
163: {
164:   PetscFunctionBegin;
166:   PetscAssertPointer(snes, 2);
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, this MUST have an attached `DM`/`PetscDS`, so that the exact solution can be computed

180:   Level: intermediate

182: .seealso: `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
183: @*/
184: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
185: {
186:   PetscFunctionBegin;
189:   adaptor->snes = snes;
190:   PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
191:   PetscFunctionReturn(PETSC_SUCCESS);
192: }

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

197:   Not Collective

199:   Input Parameter:
200: . adaptor - The `DMAdaptor` object

202:   Output Parameter:
203: . num - The number of adaptations

205:   Level: intermediate

207: .seealso: `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
208: @*/
209: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
210: {
211:   PetscFunctionBegin;
213:   PetscAssertPointer(num, 2);
214:   *num = adaptor->numSeq;
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: /*@
219:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

221:   Not Collective

223:   Input Parameters:
224: + adaptor - The `DMAdaptor` object
225: - num     - The number of adaptations

227:   Level: intermediate

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

239: /*@
240:   DMAdaptorSetUp - After the solver is specified, creates data structures for controlling adaptivity

242:   Collective

244:   Input Parameter:
245: . adaptor - The `DMAdaptor` object

247:   Level: beginner

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

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

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

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

284: static PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
285: {
286:   DM           plex;
287:   PetscDS      prob;
288:   PetscObject  obj;
289:   PetscClassId id;
290:   PetscBool    isForest;

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

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

360: static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
361: {
362:   PetscReal time = 0.0;
363:   Mat       interp;
364:   void     *ctx;

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

388: static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
389: {
390:   PetscDS      prob;
391:   PetscObject  obj;
392:   PetscClassId id;

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

401:     PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
402:     /* Restore original limiter */
403:     PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));

405:     PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
406:     PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
407:     PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
408:   }
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: /*
413:   DMAdaptorSimpleErrorIndicator - Use the integrated gradient as an error indicator in the `DMAdaptor`

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

424:   Output Parameter:
425: . errInd   - The error indicator

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

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

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

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

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

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

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

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

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

526: 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[])
527: {
528:   PetscInt i, j;

530:   for (i = 0; i < dim; ++i) {
531:     for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
532:   }
533: }

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

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

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

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

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

594:       PetscCall(DMConvert(dm, DMPLEX, &plex));
595:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
596:       PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));

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

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

642:       PetscCall(PetscMalloc(1, &funcs));
643:       funcs[0] = identityFunc;

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

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

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

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

772:   Not Collective

774:   Input Parameters:
775: + adaptor  - The `DMAdaptor` object
776: . x        - The global approximate solution
777: - strategy - The adaptation strategy, see `DMAdaptationStrategy`

779:   Output Parameters:
780: + adm - The adapted `DM`
781: - ax  - The adapted solution

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

789:   Level: intermediate

791: .seealso: `DMAdaptor`, `DMAdaptationStrategy`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`
792: @*/
793: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
794: {
795:   PetscFunctionBegin;
796:   switch (strategy) {
797:   case DM_ADAPTATION_INITIAL:
798:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
799:     break;
800:   case DM_ADAPTATION_SEQUENTIAL:
801:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
802:     break;
803:   default:
804:     SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
805:   }
806:   PetscFunctionReturn(PETSC_SUCCESS);
807: }