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: {

 18:   DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);
 19:   return(0);
 20: }

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

 25:   Collective

 27:   Input Parameter:
 28: . comm - The communicator for the DMAdaptor object

 30:   Output Parameter:
 31: . adaptor   - The DMAdaptor object

 33:   Level: beginner

 35: .seealso: DMAdaptorDestroy(), DMAdaptorAdapt()
 36: @*/
 37: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
 38: {
 39:   VecTaggerBox     refineBox, coarsenBox;
 40:   PetscErrorCode   ierr;

 44:   PetscSysInitializePackage();
 45:   PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);
 46:   PetscNew(&(*adaptor)->metricCtx);

 48:   (*adaptor)->monitor = PETSC_FALSE;
 49:   (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
 50:   (*adaptor)->numSeq = 1;
 51:   (*adaptor)->Nadapt = -1;
 52:   (*adaptor)->refinementFactor = 2.0;
 53:   (*adaptor)->metricCtx->h_min = 1.0e-05;
 54:   (*adaptor)->metricCtx->h_max = 1.0e+05;
 55:   (*adaptor)->metricCtx->a_max = 1.0e+05;
 56:   (*adaptor)->metricCtx->p = 1.0;
 57:   (*adaptor)->metricCtx->isotropic = PETSC_FALSE;
 58:   (*adaptor)->metricCtx->restrictAnisotropyFirst = PETSC_FALSE;
 59:   (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
 60:   refineBox.min = refineBox.max = PETSC_MAX_REAL;
 61:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);
 62:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");
 63:   VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);
 64:   VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);
 65:   coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
 66:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);
 67:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");
 68:   VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);
 69:   VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);
 70:   return(0);
 71: }

 73: /*@
 74:   DMAdaptorDestroy - Destroys a DMAdaptor object

 76:   Collective on DMAdaptor

 78:   Input Parameter:
 79: . adaptor - The DMAdaptor object

 81:   Level: beginner

 83: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
 84: @*/
 85: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
 86: {

 90:   if (!*adaptor) return(0);
 92:   if (--((PetscObject)(*adaptor))->refct > 0) {
 93:     *adaptor = NULL;
 94:     return(0);
 95:   }
 96:   PetscFree((*adaptor)->metricCtx);
 97:   VecTaggerDestroy(&(*adaptor)->refineTag);
 98:   VecTaggerDestroy(&(*adaptor)->coarsenTag);
 99:   PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);
100:   PetscHeaderDestroy(adaptor);
101:   return(0);
102: }

104: /*@
105:   DMAdaptorSetFromOptions - Sets a DMAdaptor object from options

107:   Collective on DMAdaptor

109:   Input Parameter:
110: . adaptor - The DMAdaptor object

112:   Options Database Keys:
113: + -adaptor_monitor <bool>        : Monitor the adaptation process
114: . -adaptor_sequence_num <num>    : Number of adaptations to generate an optimal grid
115: . -adaptor_target_num <num>      : Set the target number of vertices N_adapt, -1 for automatic determination
116: . -adaptor_refinement_factor <r> : Set r such that N_adapt = r^dim N_orig
117: . -adaptor_metric_h_min <min>    : Set the minimum tolerated metric magnitude
118: . -adaptor_metric_h_max <max>    : Set the maximum tolerated metric magnitude
119: . -adaptor_metric_a_max <max>    : Set the maximum tolerated anisotropy
120: - -adaptor_metric_p <p>          : Set the L-p normalization order

122:   Level: beginner

124: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
125: @*/
126: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
127: {

131:   PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");
132:   PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
133:   PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
134:   PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
135:   PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
136:   PetscOptionsReal("-adaptor_metric_h_min", "Set the minimum eigenvalue of Hessian (sqr max edge length)", "DMAdaptor", adaptor->metricCtx->h_min, &adaptor->metricCtx->h_min, NULL);
137:   PetscOptionsReal("-adaptor_metric_h_max", "Set the maximum eigenvalue of Hessian (sqr min edge length)", "DMAdaptor", adaptor->metricCtx->h_max, &adaptor->metricCtx->h_max, NULL);
138:   PetscOptionsReal("-adaptor_metric_a_max", "Set the maximum tolerated anisotropy", "DMAdaptor", adaptor->metricCtx->a_max, &adaptor->metricCtx->a_max, NULL);
139:   PetscOptionsReal("-adaptor_metric_p", "Set the metric L-p normalization order", "DMAdaptor", adaptor->metricCtx->p, &adaptor->metricCtx->p, NULL);
140:   PetscOptionsEnd();
141:   VecTaggerSetFromOptions(adaptor->refineTag);
142:   VecTaggerSetFromOptions(adaptor->coarsenTag);
143:   return(0);
144: }

146: /*@
147:   DMAdaptorView - Views a DMAdaptor object

149:   Collective on DMAdaptor

151:   Input Parameters:
152: + adaptor     - The DMAdaptor object
153: - viewer - The PetscViewer object

155:   Level: beginner

157: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
158: @*/
159: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
160: {

164:   PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
165:   PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
166:   PetscViewerASCIIPrintf(viewer, "  sequence length: %D\n", adaptor->numSeq);
167:   VecTaggerView(adaptor->refineTag,  viewer);
168:   VecTaggerView(adaptor->coarsenTag, viewer);
169:   return(0);
170: }

172: /*@
173:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

175:   Not collective

177:   Input Parameter:
178: . adaptor   - The DMAdaptor object

180:   Output Parameter:
181: . snes - The solver

183:   Level: intermediate

185: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
186: @*/
187: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
188: {
192:   *snes = adaptor->snes;
193:   return(0);
194: }

196: /*@
197:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

199:   Not collective

201:   Input Parameters:
202: + adaptor   - The DMAdaptor object
203: - snes - The solver

205:   Level: intermediate

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

209: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
210: @*/
211: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
212: {

218:   adaptor->snes = snes;
219:   SNESGetDM(adaptor->snes, &adaptor->idm);
220:   return(0);
221: }

223: /*@
224:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations

226:   Not collective

228:   Input Parameter:
229: . adaptor - The DMAdaptor object

231:   Output Parameter:
232: . num - The number of adaptations

234:   Level: intermediate

236: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
237: @*/
238: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
239: {
243:   *num = adaptor->numSeq;
244:   return(0);
245: }

247: /*@
248:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

250:   Not collective

252:   Input Parameters:
253: + adaptor - The DMAdaptor object
254: - num - The number of adaptations

256:   Level: intermediate

258: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
259: @*/
260: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
261: {
264:   adaptor->numSeq = num;
265:   return(0);
266: }

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

271:   Collective on DMAdaptor

273:   Input Parameters:
274: . adaptor - The DMAdaptor object

276:   Level: beginner

278: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
279: @*/
280: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
281: {
282:   PetscDS        prob;
283:   PetscInt       Nf, f;

287:   DMGetDS(adaptor->idm, &prob);
288:   VecTaggerSetUp(adaptor->refineTag);
289:   VecTaggerSetUp(adaptor->coarsenTag);
290:   PetscDSGetNumFields(prob, &Nf);
291:   PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
292:   for (f = 0; f < Nf; ++f) {
293:     PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
294:     /* TODO Have a flag that forces projection rather than using the exact solution */
295:     if (adaptor->exactSol[0]) {DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);}
296:   }
297:   return(0);
298: }

300: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
301: {
303:   *tfunc = adaptor->ops->transfersolution;
304:   return(0);
305: }

307: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
308: {
310:   adaptor->ops->transfersolution = tfunc;
311:   return(0);
312: }

314: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
315: {
316:   DM             plex;
317:   PetscDS        prob;
318:   PetscObject    obj;
319:   PetscClassId   id;
320:   PetscBool      isForest;

324:   DMConvert(adaptor->idm, DMPLEX, &plex);
325:   DMGetDS(adaptor->idm, &prob);
326:   PetscDSGetDiscretization(prob, 0, &obj);
327:   PetscObjectGetClassId(obj, &id);
328:   DMIsForest(adaptor->idm, &isForest);
329:   if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
330:     if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
331: #if defined(PETSC_HAVE_PRAGMATIC)
332:     else          {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
333: #else
334:     else          {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
335: #endif
336:   }
337:   if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
338:   else                       {adaptor->femType = PETSC_TRUE;}
339:   if (adaptor->femType) {
340:     /* Compute local solution bc */
341:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
342:   } else {
343:     PetscFV      fvm = (PetscFV) obj;
344:     PetscLimiter noneLimiter;
345:     Vec          grad;

347:     PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
348:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
349:     /* Use no limiting when reconstructing gradients for adaptivity */
350:     PetscFVGetLimiter(fvm, &adaptor->limiter);
351:     PetscObjectReference((PetscObject) adaptor->limiter);
352:     PetscLimiterCreate(PetscObjectComm((PetscObject) fvm), &noneLimiter);
353:     PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
354:     PetscFVSetLimiter(fvm, noneLimiter);
355:     /* Get FVM data */
356:     DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
357:     VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
358:     VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
359:     /* Compute local solution bc */
360:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
361:     /* Compute gradients */
362:     DMCreateGlobalVector(adaptor->gradDM, &grad);
363:     DMPlexReconstructGradientsFVM(plex, locX, grad);
364:     DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
365:     DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
366:     DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
367:     VecDestroy(&grad);
368:     VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
369:   }
370:   DMDestroy(&plex);
371:   return(0);
372: }

374: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
375: {
376:   PetscReal      time = 0.0;
377:   Mat            interp;
378:   void          *ctx;

382:   DMGetApplicationContext(dm, &ctx);
383:   if (adaptor->ops->transfersolution) {
384:     (*adaptor->ops->transfersolution)(adaptor, dm, x, adm, ax, ctx);
385:   } else {
386:     switch (adaptor->adaptCriterion) {
387:     case DM_ADAPTATION_LABEL:
388:       DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
389:       break;
390:     case DM_ADAPTATION_REFINE:
391:     case DM_ADAPTATION_METRIC:
392:       DMCreateInterpolation(dm, adm, &interp, NULL);
393:       MatInterpolate(interp, x, ax);
394:       DMInterpolate(dm, interp, adm);
395:       MatDestroy(&interp);
396:       break;
397:     default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %D", adaptor->adaptCriterion);
398:     }
399:   }
400:   return(0);
401: }

403: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
404: {
405:   PetscDS        prob;
406:   PetscObject    obj;
407:   PetscClassId   id;

411:   DMGetDS(adaptor->idm, &prob);
412:   PetscDSGetDiscretization(prob, 0, &obj);
413:   PetscObjectGetClassId(obj, &id);
414:   if (id == PETSCFV_CLASSID) {
415:     PetscFV fvm = (PetscFV) obj;

417:     PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
418:     /* Restore original limiter */
419:     PetscFVSetLimiter(fvm, adaptor->limiter);

421:     VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
422:     VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
423:     DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
424:   }
425:   return(0);
426: }

428: /*
429:   DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator

431:   Input Parameters:
432: + adaptor  - The DMAdaptor object
433: . dim      - The topological dimension
434: . cell     - The cell
435: . field    - The field integrated over the cell
436: . gradient - The gradient integrated over the cell
437: . cg       - A PetscFVCellGeom struct
438: - ctx      - A user context

440:   Output Parameter:
441: . errInd   - The error indicator

443: .seealso: DMAdaptorComputeErrorIndicator()
444: */
445: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
446: {
447:   PetscReal err = 0.;
448:   PetscInt  c, d;

451:   for (c = 0; c < Nc; c++) {
452:     for (d = 0; d < dim; ++d) {
453:       err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
454:     }
455:   }
456:   *errInd = cg->volume * err;
457:   return(0);
458: }

460: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
461: {
462:   PetscDS         prob;
463:   PetscObject     obj;
464:   PetscClassId    id;
465:   void           *ctx;
466:   PetscQuadrature quad;
467:   PetscInt        dim, d, cdim, Nc;
468:   PetscErrorCode  ierr;

471:   *errInd = 0.;
472:   DMGetDimension(plex, &dim);
473:   DMGetCoordinateDim(plex, &cdim);
474:   DMGetApplicationContext(plex, &ctx);
475:   DMGetDS(plex, &prob);
476:   PetscDSGetDiscretization(prob, 0, &obj);
477:   PetscObjectGetClassId(obj, &id);
478:   if (id == PETSCFV_CLASSID) {
479:     const PetscScalar *pointSols;
480:     const PetscScalar *pointSol;
481:     const PetscScalar *pointGrad;
482:     PetscFVCellGeom   *cg;

484:     PetscFVGetNumComponents((PetscFV) obj, &Nc);
485:     VecGetArrayRead(locX, &pointSols);
486:     DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
487:     DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
488:     DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
489:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
490:     VecRestoreArrayRead(locX, &pointSols);
491:   } else {
492:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
493:     PetscFVCellGeom  cg;
494:     PetscFEGeom      fegeom;
495:     const PetscReal *quadWeights;
496:     PetscReal       *coords;
497:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

499:     fegeom.dim      = dim;
500:     fegeom.dimEmbed = cdim;
501:     PetscDSGetNumFields(prob, &Nf);
502:     PetscFEGetQuadrature((PetscFE) obj, &quad);
503:     DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
504:     PetscFEGetDimension((PetscFE) obj, &Nb);
505:     PetscFEGetNumComponents((PetscFE) obj, &Nc);
506:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
507:     PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);
508:     PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
509:     DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
510:     DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
511:     PetscArrayzero(gradient, cdim*Nc);
512:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
513:       PetscInt qc = 0, q;

515:       PetscDSGetDiscretization(prob, f, &obj);
516:       PetscArrayzero(interpolant,Nc);
517:       PetscArrayzero(interpolantGrad, cdim*Nc);
518:       for (q = 0; q < Nq; ++q) {
519:         PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, 1, x, &fegeom, q, interpolant, interpolantGrad);
520:         for (fc = 0; fc < Nc; ++fc) {
521:           const PetscReal wt = quadWeights[q*qNc+qc+fc];

523:           field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
524:           for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
525:         }
526:       }
527:       fieldOffset += Nb;
528:       qc          += Nc;
529:     }
530:     PetscFree2(interpolant, interpolantGrad);
531:     DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
532:     for (fc = 0; fc < Nc; ++fc) {
533:       field[fc] /= cg.volume;
534:       for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
535:     }
536:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);
537:     PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
538:   }
539:   return(0);
540: }

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

552:   DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
553:   VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
554:   PetscObjectGetComm((PetscObject) adaptor, &comm);
555:   DMGetDimension(adaptor->idm, &dim);
556:   DMGetCoordinateDim(adaptor->idm, &coordDim);
557:   DMGetApplicationContext(adaptor->idm, &ctx);
558:   DMGetDS(adaptor->idm, &prob);
559:   PetscDSGetNumFields(prob, &numFields);
560:   if (numFields == 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fields is zero!");

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

570:     DMGetLocalVector(dm, &locX);
571:     DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
572:     DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
573:     DMAdaptorPreAdapt(adaptor, locX);
574:     if (doSolve) {
575:       SNES snes;

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

603:       DMConvert(dm, DMPLEX, &plex);
604:       DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
605:       DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);

607:       VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
608:       VecSetUp(errVec);
609:       VecGetArray(errVec, &errArray);
610:       VecGetArrayRead(locX, &pointSols);
611:       for (c = cStart; c < cEnd; ++c) {
612:         PetscReal errInd;

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

651:       /*     Compute vertexwise gradients from cellwise gradients */
652:       DMClone(dm, &dmGrad);
653:       DMClone(dm, &dmHess);
654:       for (f = 0; f < numFields; ++f) {
655:         PetscFE         fe, feGrad, feHess;
656:         PetscDualSpace  Q;
657:         DM              K;
658:         PetscQuadrature q;
659:         PetscInt        Nc, qorder;
660:         const char     *prefix;

662:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
663:         PetscFEGetNumComponents(fe, &Nc);
664:         PetscFEGetDualSpace(fe, &Q);
665:         PetscDualSpaceGetDM(Q, &K);
666:         DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
667:         PetscFEGetQuadrature(fe, &q);
668:         PetscQuadratureGetOrder(q, &qorder);
669:         PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);
670:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feGrad);
671:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feHess);
672:         DMSetField(dmGrad, f, NULL, (PetscObject)feGrad);
673:         DMSetField(dmHess, f, NULL, (PetscObject)feHess);
674:         DMCreateDS(dmGrad);
675:         DMCreateDS(dmHess);
676:         PetscFEDestroy(&feGrad);
677:         PetscFEDestroy(&feHess);
678:       }
679:       DMGetGlobalVector(dmGrad, &xGrad);
680:       VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view");
681:       DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
682:       VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
683:       /*     Compute vertexwise Hessians from cellwise Hessians */
684:       DMGetGlobalVector(dmHess, &xHess);
685:       DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
686:       VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
687:       DMRestoreGlobalVector(dmGrad, &xGrad);
688:       DMDestroy(&dmGrad);
689:       /*     Set target metric complexity */
690:       N    = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
691:       if (adaptor->monitor) {PetscPrintf(PETSC_COMM_SELF, "N_orig: %D N_adapt: %g\n", vEnd - vStart, N);}
692:       adaptor->metricCtx->targetComplexity = (PetscReal) N;
693:       DMClone(dm, &dmMetric);
694:       DMSetApplicationContext(dmMetric, adaptor->metricCtx);
695:       /*     Compute L-p normalized metric */
696:       DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, &metric);
697:       DMRestoreGlobalVector(dmHess, &xHess);
698:       DMDestroy(&dmHess);
699:       /*     Adapt DM from metric */
700:       DMGetLabel(dm, "marker", &bdLabel);
701:       DMAdaptMetric(dm, metric, bdLabel, &odm);
702:       adapted = PETSC_TRUE;
703:       /*     Cleanup */
704:       DMRestoreLocalVector(dmMetric, &metric);
705:       DMDestroy(&dmMetric);
706:     }
707:     break;
708:     default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
709:     }
710:     DMAdaptorPostAdapt(adaptor);
711:     DMRestoreLocalVector(dm, &locX);
712:     /* If DM was adapted, replace objects and recreate solution */
713:     if (adapted) {
714:       const char *name;

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

753: /*@
754:   DMAdaptorAdapt - Creates a new DM that is adapted to the problem

756:   Not collective

758:   Input Parameters:
759: + adaptor  - The DMAdaptor object
760: . x        - The global approximate solution
761: - strategy - The adaptation strategy

763:   Output Parameters:
764: + adm - The adapted DM
765: - ax  - The adapted solution

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

773:   Note: The available adaptation strategies are:
774: $ 1) Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
775: $ 2) Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
776: $ 3) Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid

778:   Level: intermediate

780: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
781: @*/
782: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
783: {

787:   switch (strategy)
788:   {
789:   case DM_ADAPTATION_INITIAL:
790:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
791:     break;
792:   case DM_ADAPTATION_SEQUENTIAL:
793:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
794:     break;
795:   default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
796:   }
797:   return(0);
798: }