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: {
 16:   DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);
 17:   return 0;
 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;

 40:   PetscSysInitializePackage();
 41:   PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);

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

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

 65:   Collective

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

 70:   Level: beginner

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

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

 92:   Collective

 94:   Input Parameter:
 95: . adaptor - The `DMAdaptor` object

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

103:   Level: beginner

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

120: /*@
121:    DMAdaptorView - Views a `DMAdaptor` object

123:    Collective

125:    Input Parameters:
126: +  adaptor - The `DMAdaptor` object
127: -  viewer - The `PetscViewer` object

129:   Level: beginner

131: .seealso: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
132: @*/
133: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
134: {
135:   PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer);
136:   PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
137:   PetscViewerASCIIPrintf(viewer, "  sequence length: %" PetscInt_FMT "\n", adaptor->numSeq);
138:   VecTaggerView(adaptor->refineTag, viewer);
139:   VecTaggerView(adaptor->coarsenTag, viewer);
140:   return 0;
141: }

143: /*@
144:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

146:   Not collective

148:   Input Parameter:
149: . adaptor   - The `DMAdaptor` object

151:   Output Parameter:
152: . snes - The solver

154:   Level: intermediate

156: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
157: @*/
158: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
159: {
162:   *snes = adaptor->snes;
163:   return 0;
164: }

166: /*@
167:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

169:   Not collective

171:   Input Parameters:
172: + adaptor   - The `DMAdaptor` object
173: - snes - The solver

175:   Level: intermediate

177:   Note:
178:   The solver MUST have an attached `DM`/`DS`, so that we know the exact solution

180: .seealso: `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
181: @*/
182: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
183: {
186:   adaptor->snes = snes;
187:   SNESGetDM(adaptor->snes, &adaptor->idm);
188:   return 0;
189: }

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

194:   Not collective

196:   Input Parameter:
197: . adaptor - The `DMAdaptor` object

199:   Output Parameter:
200: . num - The number of adaptations

202:   Level: intermediate

204: .seealso: `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
205: @*/
206: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
207: {
210:   *num = adaptor->numSeq;
211:   return 0;
212: }

214: /*@
215:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

217:   Not collective

219:   Input Parameters:
220: + adaptor - The `DMAdaptor` object
221: - num - The number of adaptations

223:   Level: intermediate

225: .seealso: `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
226: @*/
227: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
228: {
230:   adaptor->numSeq = num;
231:   return 0;
232: }

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

237:   Collective

239:   Input Parameters:
240: . adaptor - The `DMAdaptor` object

242:   Level: beginner

244: .seealso: `DMAdaptorCreate()`, `DMAdaptorAdapt()`
245: @*/
246: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
247: {
248:   PetscDS  prob;
249:   PetscInt Nf, f;

251:   DMGetDS(adaptor->idm, &prob);
252:   VecTaggerSetUp(adaptor->refineTag);
253:   VecTaggerSetUp(adaptor->coarsenTag);
254:   PetscDSGetNumFields(prob, &Nf);
255:   PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
256:   for (f = 0; f < Nf; ++f) {
257:     PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
258:     /* TODO Have a flag that forces projection rather than using the exact solution */
259:     if (adaptor->exactSol[0]) DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);
260:   }
261:   return 0;
262: }

264: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
265: {
266:   *tfunc = adaptor->ops->transfersolution;
267:   return 0;
268: }

270: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
271: {
272:   adaptor->ops->transfersolution = tfunc;
273:   return 0;
274: }

276: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
277: {
278:   DM           plex;
279:   PetscDS      prob;
280:   PetscObject  obj;
281:   PetscClassId id;
282:   PetscBool    isForest;

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

324:     PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
325:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
326:     /* Use no limiting when reconstructing gradients for adaptivity */
327:     PetscFVGetLimiter(fvm, &adaptor->limiter);
328:     PetscObjectReference((PetscObject)adaptor->limiter);
329:     PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter);
330:     PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
331:     PetscFVSetLimiter(fvm, noneLimiter);
332:     /* Get FVM data */
333:     DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
334:     VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
335:     VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
336:     /* Compute local solution bc */
337:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
338:     /* Compute gradients */
339:     DMCreateGlobalVector(adaptor->gradDM, &grad);
340:     DMPlexReconstructGradientsFVM(plex, locX, grad);
341:     DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
342:     DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
343:     DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
344:     VecDestroy(&grad);
345:     VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
346:   }
347:   DMDestroy(&plex);
348:   return 0;
349: }

351: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
352: {
353:   PetscReal time = 0.0;
354:   Mat       interp;
355:   void     *ctx;

357:   DMGetApplicationContext(dm, &ctx);
358:   if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
359:   else {
360:     switch (adaptor->adaptCriterion) {
361:     case DM_ADAPTATION_LABEL:
362:       DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
363:       break;
364:     case DM_ADAPTATION_REFINE:
365:     case DM_ADAPTATION_METRIC:
366:       DMCreateInterpolation(dm, adm, &interp, NULL);
367:       MatInterpolate(interp, x, ax);
368:       DMInterpolate(dm, interp, adm);
369:       MatDestroy(&interp);
370:       break;
371:     default:
372:       SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
373:     }
374:   }
375:   return 0;
376: }

378: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
379: {
380:   PetscDS      prob;
381:   PetscObject  obj;
382:   PetscClassId id;

384:   DMGetDS(adaptor->idm, &prob);
385:   PetscDSGetDiscretization(prob, 0, &obj);
386:   PetscObjectGetClassId(obj, &id);
387:   if (id == PETSCFV_CLASSID) {
388:     PetscFV fvm = (PetscFV)obj;

390:     PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
391:     /* Restore original limiter */
392:     PetscFVSetLimiter(fvm, adaptor->limiter);

394:     VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
395:     VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
396:     DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
397:   }
398:   return 0;
399: }

401: /*
402:   DMAdaptorSimpleErrorIndicator - Use the integrated gradient as an error indicator

404:   Input Parameters:
405: + adaptor  - The `DMAdaptor` object
406: . dim      - The topological dimension
407: . cell     - The cell
408: . field    - The field integrated over the cell
409: . gradient - The gradient integrated over the cell
410: . cg       - A `PetscFVCellGeom` struct
411: - ctx      - A user context

413:   Output Parameter:
414: . errInd   - The error indicator

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

419: .seealso: `DMAdaptor`, `DMAdaptorComputeErrorIndicator()`
420: */
421: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
422: {
423:   PetscReal err = 0.;
424:   PetscInt  c, d;

427:   for (c = 0; c < Nc; c++) {
428:     for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
429:   }
430:   *errInd = cg->volume * err;
431:   return 0;
432: }

434: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
435: {
436:   PetscDS         prob;
437:   PetscObject     obj;
438:   PetscClassId    id;
439:   void           *ctx;
440:   PetscQuadrature quad;
441:   PetscInt        dim, d, cdim, Nc;

443:   *errInd = 0.;
444:   DMGetDimension(plex, &dim);
445:   DMGetCoordinateDim(plex, &cdim);
446:   DMGetApplicationContext(plex, &ctx);
447:   DMGetDS(plex, &prob);
448:   PetscDSGetDiscretization(prob, 0, &obj);
449:   PetscObjectGetClassId(obj, &id);
450:   if (id == PETSCFV_CLASSID) {
451:     const PetscScalar *pointSols;
452:     const PetscScalar *pointSol;
453:     const PetscScalar *pointGrad;
454:     PetscFVCellGeom   *cg;

456:     PetscFVGetNumComponents((PetscFV)obj, &Nc);
457:     VecGetArrayRead(locX, &pointSols);
458:     DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol);
459:     DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad);
460:     DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
461:     PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
462:     VecRestoreArrayRead(locX, &pointSols);
463:   } else {
464:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
465:     PetscFVCellGeom  cg;
466:     PetscFEGeom      fegeom;
467:     const PetscReal *quadWeights;
468:     PetscReal       *coords;
469:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

471:     fegeom.dim      = dim;
472:     fegeom.dimEmbed = cdim;
473:     PetscDSGetNumFields(prob, &Nf);
474:     PetscFEGetQuadrature((PetscFE)obj, &quad);
475:     DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
476:     PetscFEGetDimension((PetscFE)obj, &Nb);
477:     PetscFEGetNumComponents((PetscFE)obj, &Nc);
478:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
479:     PetscMalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ);
480:     PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad);
481:     DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
482:     DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
483:     PetscArrayzero(gradient, cdim * Nc);
484:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
485:       PetscInt qc = 0, q;

487:       PetscDSGetDiscretization(prob, f, &obj);
488:       PetscArrayzero(interpolant, Nc);
489:       PetscArrayzero(interpolantGrad, cdim * Nc);
490:       for (q = 0; q < Nq; ++q) {
491:         PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad);
492:         for (fc = 0; fc < Nc; ++fc) {
493:           const PetscReal wt = quadWeights[q * qNc + qc + fc];

495:           field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
496:           for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
497:         }
498:       }
499:       fieldOffset += Nb;
500:       qc += Nc;
501:     }
502:     PetscFree2(interpolant, interpolantGrad);
503:     DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
504:     for (fc = 0; fc < Nc; ++fc) {
505:       field[fc] /= cg.volume;
506:       for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
507:     }
508:     PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, field, gradient, &cg, errInd, ctx);
509:     PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ);
510:   }
511:   return 0;
512: }

514: 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[])
515: {
516:   PetscInt i, j;

518:   for (i = 0; i < dim; ++i) {
519:     for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
520:   }
521: }

523: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
524: {
525:   PetscDS  prob;
526:   void    *ctx;
527:   MPI_Comm comm;
528:   PetscInt numAdapt = adaptor->numSeq, adaptIter;
529:   PetscInt dim, coordDim, numFields, cStart, cEnd, c;

531:   DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
532:   VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
533:   PetscObjectGetComm((PetscObject)adaptor, &comm);
534:   DMGetDimension(adaptor->idm, &dim);
535:   DMGetCoordinateDim(adaptor->idm, &coordDim);
536:   DMGetApplicationContext(adaptor->idm, &ctx);
537:   DMGetDS(adaptor->idm, &prob);
538:   PetscDSGetNumFields(prob, &numFields);

541:   /* Adapt until nothing changes */
542:   /* Adapt for a specified number of iterates */
543:   for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));
544:   for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
545:     PetscBool adapted = PETSC_FALSE;
546:     DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
547:     Vec       x       = adaptIter ? *ax : inx, locX, ox;

549:     DMGetLocalVector(dm, &locX);
550:     DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
551:     DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
552:     DMAdaptorPreAdapt(adaptor, locX);
553:     if (doSolve) {
554:       SNES snes;

556:       DMAdaptorGetSolver(adaptor, &snes);
557:       SNESSolve(snes, NULL, adaptIter ? *ax : x);
558:     }
559:     /* DMAdaptorMonitor(adaptor);
560:        Print iterate, memory used, DM, solution */
561:     switch (adaptor->adaptCriterion) {
562:     case DM_ADAPTATION_REFINE:
563:       DMRefine(dm, comm, &odm);
565:       adapted = PETSC_TRUE;
566:       break;
567:     case DM_ADAPTATION_LABEL: {
568:       /* Adapt DM
569:            Create local solution
570:            Reconstruct gradients (FVM) or solve adjoint equation (FEM)
571:            Produce cellwise error indicator */
572:       DM                 plex;
573:       DMLabel            adaptLabel;
574:       IS                 refineIS, coarsenIS;
575:       Vec                errVec;
576:       PetscScalar       *errArray;
577:       const PetscScalar *pointSols;
578:       PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
579:       PetscInt           nRefine, nCoarsen;

581:       DMConvert(dm, DMPLEX, &plex);
582:       DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
583:       DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);

585:       VecCreateMPI(PetscObjectComm((PetscObject)adaptor), cEnd - cStart, PETSC_DETERMINE, &errVec);
586:       VecSetUp(errVec);
587:       VecGetArray(errVec, &errArray);
588:       VecGetArrayRead(locX, &pointSols);
589:       for (c = cStart; c < cEnd; ++c) {
590:         PetscReal errInd;

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

629:       PetscMalloc(1, &funcs);
630:       funcs[0] = identityFunc;

632:       /*     Setup finite element spaces */
633:       DMClone(dm, &dmGrad);
634:       DMClone(dm, &dmHess);
636:       for (f = 0; f < numFields; ++f) {
637:         PetscFE         fe, feGrad, feHess;
638:         PetscDualSpace  Q;
639:         PetscSpace      space;
640:         DM              K;
641:         PetscQuadrature q;
642:         PetscInt        Nc, qorder, p;
643:         const char     *prefix;

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

719:       PetscObjectGetName((PetscObject)dm, &name);
720:       PetscObjectSetName((PetscObject)odm, name);
721:       /* Reconfigure solver */
722:       SNESReset(adaptor->snes);
723:       SNESSetDM(adaptor->snes, odm);
724:       DMAdaptorSetSolver(adaptor, adaptor->snes);
725:       DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
726:       SNESSetFromOptions(adaptor->snes);
727:       /* Transfer system */
728:       DMCopyDisc(dm, odm);
729:       /* Transfer solution */
730:       DMCreateGlobalVector(odm, &ox);
731:       PetscObjectGetName((PetscObject)x, &name);
732:       PetscObjectSetName((PetscObject)ox, name);
733:       DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
734:       /* Cleanup adaptivity info */
735:       if (adaptIter > 0) PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));
736:       DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
737:       DMDestroy(&dm);
738:       VecDestroy(&x);
739:       *adm = odm;
740:       *ax  = ox;
741:     } else {
742:       *adm      = dm;
743:       *ax       = x;
744:       adaptIter = numAdapt;
745:     }
746:     if (adaptIter < numAdapt - 1) {
747:       DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
748:       VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
749:     }
750:   }
751:   DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
752:   VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
753:   return 0;
754: }

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

759:   Not collective

761:   Input Parameters:
762: + adaptor  - The `DMAdaptor` object
763: . x        - The global approximate solution
764: - strategy - The adaptation strategy

766:   Output Parameters:
767: + adm - The adapted `DM`
768: - ax  - The adapted solution

770:   Options database Keys:
771: + -snes_adapt <strategy> - initial, sequential, multigrid
772: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
773: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
774: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement

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

781:   Level: intermediate

783: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
784: @*/
785: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
786: {
787:   switch (strategy) {
788:   case DM_ADAPTATION_INITIAL:
789:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
790:     break;
791:   case DM_ADAPTATION_SEQUENTIAL:
792:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
793:     break;
794:   default:
795:     SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
796:   }
797:   return 0;
798: }