Actual source code: dmadapt.c

petsc-master 2019-11-13
Report Typos and Errors
  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;

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

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

 68: /*@
 69:   DMAdaptorDestroy - Destroys a DMAdaptor object

 71:   Collective on DMAdaptor

 73:   Input Parameter:
 74: . adaptor - The DMAdaptor object

 76:   Level: beginner

 78: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
 79: @*/
 80: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
 81: {

 85:   if (!*adaptor) return(0);
 87:   if (--((PetscObject)(*adaptor))->refct > 0) {
 88:     *adaptor = NULL;
 89:     return(0);
 90:   }
 91:   VecTaggerDestroy(&(*adaptor)->refineTag);
 92:   VecTaggerDestroy(&(*adaptor)->coarsenTag);
 93:   PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);
 94:   PetscHeaderDestroy(adaptor);
 95:   return(0);
 96: }

 98: /*@
 99:   DMAdaptorSetFromOptions - Sets a DMAdaptor object from options

101:   Collective on DMAdaptor

103:   Input Parameters:
104: . adaptor - The DMAdaptor object

106:   Options Database Keys:
107: + -adaptor_monitor <bool>        : Monitor the adaptation process
108: . -adaptor_sequence_num <num>    : Number of adaptations to generate an optimal grid
109: . -adaptor_target_num <num>      : Set the target number of vertices N_adapt, -1 for automatic determination
110: . -adaptor_refinement_factor <r> : Set r such that N_adapt = r^dim N_orig
111: . -adaptor_metric_h_min <min>    : Set the minimum eigenvalue of Hessian (sqr max edge length)
112: - -adaptor_metric_h_max <max>    : Set the maximum eigenvalue of Hessian (sqr min edge length)

114:   Level: beginner

116: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
117: @*/
118: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
119: {

123:   PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");
124:   PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
125:   PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
126:   PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
127:   PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
128:   PetscOptionsReal("-adaptor_metric_h_min", "Set the minimum eigenvalue of Hessian (sqr max edge length)", "DMAdaptor", adaptor->h_min, &adaptor->h_min, NULL);
129:   PetscOptionsReal("-adaptor_metric_h_max", "Set the maximum eigenvalue of Hessian (sqr min edge length)", "DMAdaptor", adaptor->h_max, &adaptor->h_max, NULL);
130:   PetscOptionsEnd();
131:   VecTaggerSetFromOptions(adaptor->refineTag);
132:   VecTaggerSetFromOptions(adaptor->coarsenTag);
133:   return(0);
134: }

136: /*@
137:   DMAdaptorView - Views a DMAdaptor object

139:   Collective on DMAdaptor

141:   Input Parameters:
142: + adaptor     - The DMAdaptor object
143: - viewer - The PetscViewer object

145:   Level: beginner

147: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
148: @*/
149: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
150: {

154:   PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
155:   PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
156:   PetscViewerASCIIPrintf(viewer, "  sequence length: %D\n", adaptor->numSeq);
157:   VecTaggerView(adaptor->refineTag,  viewer);
158:   VecTaggerView(adaptor->coarsenTag, viewer);
159:   return(0);
160: }

162: /*@
163:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

165:   Not collective

167:   Input Parameter:
168: . adaptor   - The DMAdaptor object

170:   Output Parameter:
171: . snes - The solver

173:   Level: intermediate

175: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
176: @*/
177: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
178: {
182:   *snes = adaptor->snes;
183:   return(0);
184: }

186: /*@
187:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

189:   Not collective

191:   Input Parameters:
192: + adaptor   - The DMAdaptor object
193: - snes - The solver

195:   Level: intermediate

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

199: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
200: @*/
201: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
202: {

208:   adaptor->snes = snes;
209:   SNESGetDM(adaptor->snes, &adaptor->idm);
210:   return(0);
211: }

213: /*@
214:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations

216:   Not collective

218:   Input Parameter:
219: . adaptor - The DMAdaptor object

221:   Output Parameter:
222: . num - The number of adaptations

224:   Level: intermediate

226: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
227: @*/
228: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
229: {
233:   *num = adaptor->numSeq;
234:   return(0);
235: }

237: /*@
238:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

240:   Not collective

242:   Input Parameters:
243: + adaptor - The DMAdaptor object
244: - num - The number of adaptations

246:   Level: intermediate

248: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
249: @*/
250: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
251: {
254:   adaptor->numSeq = num;
255:   return(0);
256: }

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

261:   Collective on DMAdaptor

263:   Input Parameters:
264: . adaptor - The DMAdaptor object

266:   Level: beginner

268: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
269: @*/
270: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
271: {
272:   PetscDS        prob;
273:   PetscInt       Nf, f;

277:   DMGetDS(adaptor->idm, &prob);
278:   VecTaggerSetUp(adaptor->refineTag);
279:   VecTaggerSetUp(adaptor->coarsenTag);
280:   PetscDSGetNumFields(prob, &Nf);
281:   PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
282:   for (f = 0; f < Nf; ++f) {
283:     PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
284:     /* TODO Have a flag that forces projection rather than using the exact solution */
285:     if (adaptor->exactSol[0]) {DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);}
286:   }
287:   return(0);
288: }

290: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
291: {
293:   *tfunc = adaptor->ops->transfersolution;
294:   return(0);
295: }

297: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
298: {
300:   adaptor->ops->transfersolution = tfunc;
301:   return(0);
302: }

304: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
305: {
306:   DM             plex;
307:   PetscDS        prob;
308:   PetscObject    obj;
309:   PetscClassId   id;
310:   PetscBool      isForest;

314:   DMConvert(adaptor->idm, DMPLEX, &plex);
315:   DMGetDS(adaptor->idm, &prob);
316:   PetscDSGetDiscretization(prob, 0, &obj);
317:   PetscObjectGetClassId(obj, &id);
318:   DMIsForest(adaptor->idm, &isForest);
319:   if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
320:     if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
321: #if defined(PETSC_HAVE_PRAGMATIC)
322:     else          {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
323: #else
324:     else          {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
325: #endif
326:   }
327:   if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
328:   else                       {adaptor->femType = PETSC_TRUE;}
329:   if (adaptor->femType) {
330:     /* Compute local solution bc */
331:     DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
332:   } else {
333:     PetscFV      fvm = (PetscFV) obj;
334:     PetscLimiter noneLimiter;
335:     Vec          grad;

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

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

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

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

401:   DMGetDS(adaptor->idm, &prob);
402:   PetscDSGetDiscretization(prob, 0, &obj);
403:   PetscObjectGetClassId(obj, &id);
404:   if (id == PETSCFV_CLASSID) {
405:     PetscFV fvm = (PetscFV) obj;

407:     PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
408:     /* Restore original limiter */
409:     PetscFVSetLimiter(fvm, adaptor->limiter);

411:     VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
412:     VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
413:     DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
414:   }
415:   return(0);
416: }

418: static PetscErrorCode DMAdaptorModifyHessian_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscScalar Hp[])
419: {
420:   PetscScalar   *Hpos;
421:   PetscReal     *eigs;
422:   PetscInt       i, j, k;

426:   PetscMalloc2(dim*dim, &Hpos, dim, &eigs);
427: #if 0
428:   PetscPrintf(PETSC_COMM_SELF, "H = [");
429:   for (i = 0; i < dim; ++i) {
430:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "     ");}
431:     for (j = 0; j < dim; ++j) {
432:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
433:       PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
434:     }
435:     PetscPrintf(PETSC_COMM_SELF, "\n");
436:   }
437:   PetscPrintf(PETSC_COMM_SELF, "]\n");
438: #endif
439:   /* Symmetrize */
440:   for (i = 0; i < dim; ++i) {
441:     Hpos[i*dim+i] = Hp[i*dim+i];
442:     for (j = i+1; j < dim; ++j) {
443:       Hpos[i*dim+j] = 0.5*(Hp[i*dim+j] + Hp[j*dim+i]);
444:       Hpos[j*dim+i] = Hpos[i*dim+j];
445:     }
446:   }
447: #if 0
448:   PetscPrintf(PETSC_COMM_SELF, "Hs = [");
449:   for (i = 0; i < dim; ++i) {
450:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "      ");}
451:     for (j = 0; j < dim; ++j) {
452:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
453:       PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
454:     }
455:     PetscPrintf(PETSC_COMM_SELF, "\n");
456:   }
457:   PetscPrintf(PETSC_COMM_SELF, "]\n");
458: #endif
459:   /* Compute eigendecomposition */
460:   {
461:     PetscScalar  *work;
462:     PetscBLASInt lwork;

464:     lwork = 5*dim;
465:     PetscMalloc1(5*dim, &work);
466: #if defined(PETSC_MISSING_LAPACK_GEEV)
467:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
468: #else
469:     {
470:       PetscBLASInt lierr;
471:       PetscBLASInt nb;

473:       PetscBLASIntCast(dim, &nb);
474:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
475: #if defined(PETSC_USE_COMPLEX)
476:       {
477:         PetscReal *rwork;
478:         PetscMalloc1(3*dim, &rwork);
479:         PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,rwork,&lierr));
480:         PetscFree(rwork);
481:       }
482: #else
483:       PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,&lierr));
484: #endif
485:       if (lierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
486:       PetscFPTrapPop();
487:     }
488: #endif
489:     PetscFree(work);
490:   }
491: #if 0
492:   PetscPrintf(PETSC_COMM_SELF, "L = [");
493:   for (i = 0; i < dim; ++i) {
494:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
495:     PetscPrintf(PETSC_COMM_SELF, "%g", eigs[i]);
496:   }
497:   PetscPrintf(PETSC_COMM_SELF, "]\n");
498:   PetscPrintf(PETSC_COMM_SELF, "Q = [");
499:   for (i = 0; i < dim; ++i) {
500:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "     ");}
501:     for (j = 0; j < dim; ++j) {
502:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
503:       PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
504:     }
505:     PetscPrintf(PETSC_COMM_SELF, "\n");
506:   }
507:   PetscPrintf(PETSC_COMM_SELF, "]\n");
508: #endif
509:   /* Reflect to positive orthant, enforce maximum and minimum size, \lambda \propto 1/h^2
510:        TODO get domain bounding box */
511:   for (i = 0; i < dim; ++i) eigs[i] = PetscMin(h_max, PetscMax(h_min, PetscAbsReal(eigs[i])));
512:   /* Reconstruct Hessian */
513:   for (i = 0; i < dim; ++i) {
514:     for (j = 0; j < dim; ++j) {
515:       Hp[i*dim+j] = 0.0;
516:       for (k = 0; k < dim; ++k) {
517:         Hp[i*dim+j] += Hpos[k*dim+i] * eigs[k] * Hpos[k*dim+j];
518:       }
519:     }
520:   }
521:   PetscFree2(Hpos, eigs);
522: #if 0
523:   PetscPrintf(PETSC_COMM_SELF, "H+ = [");
524:   for (i = 0; i < dim; ++i) {
525:     if (i > 0) {PetscPrintf(PETSC_COMM_SELF, "      ");}
526:     for (j = 0; j < dim; ++j) {
527:       if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
528:       PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
529:     }
530:     PetscPrintf(PETSC_COMM_SELF, "\n");
531:   }
532:   PetscPrintf(PETSC_COMM_SELF, "]\n");
533: #endif
534:   return(0);
535: }

537: static void detHFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
538:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
539:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
540:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
541: {
542:   const PetscInt p = 1;
543:   PetscReal      detH = 0.0;

545:   if      (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, u);
546:   else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, u);
547:   f0[0] = PetscPowReal(detH, p/(2.*p + dim));
548: }

550: /*
551:   DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator

553:   Input Parameters:
554: + adaptor  - The DMAdaptor object
555: . dim      - The topological dimension
556: . cell     - The cell
557: . field    - The field integrated over the cell
558: . gradient - The gradient integrated over the cell
559: . cg       - A PetscFVCellGeom struct
560: - ctx      - A user context

562:   Output Parameter:
563: . errInd   - The error indicator

565: .seealso: DMAdaptorComputeErrorIndicator()
566: */
567: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
568: {
569:   PetscReal err = 0.;
570:   PetscInt  c, d;

573:   for (c = 0; c < Nc; c++) {
574:     for (d = 0; d < dim; ++d) {
575:       err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
576:     }
577:   }
578:   *errInd = cg->volume * err;
579:   return(0);
580: }

582: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
583: {
584:   PetscDS         prob;
585:   PetscObject     obj;
586:   PetscClassId    id;
587:   void           *ctx;
588:   PetscQuadrature quad;
589:   PetscInt        dim, d, cdim, Nc;
590:   PetscErrorCode  ierr;

593:   *errInd = 0.;
594:   DMGetDimension(plex, &dim);
595:   DMGetCoordinateDim(plex, &cdim);
596:   DMGetApplicationContext(plex, &ctx);
597:   DMGetDS(plex, &prob);
598:   PetscDSGetDiscretization(prob, 0, &obj);
599:   PetscObjectGetClassId(obj, &id);
600:   if (id == PETSCFV_CLASSID) {
601:     const PetscScalar *pointSols;
602:     const PetscScalar *pointSol;
603:     const PetscScalar *pointGrad;
604:     PetscFVCellGeom   *cg;

606:     PetscFVGetNumComponents((PetscFV) obj, &Nc);
607:     VecGetArrayRead(locX, &pointSols);
608:     DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
609:     DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
610:     DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
611:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
612:     VecRestoreArrayRead(locX, &pointSols);
613:   } else {
614:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
615:     PetscFVCellGeom  cg;
616:     PetscFEGeom      fegeom;
617:     const PetscReal *quadWeights;
618:     PetscReal       *coords;
619:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

621:     fegeom.dim      = dim;
622:     fegeom.dimEmbed = cdim;
623:     PetscDSGetNumFields(prob, &Nf);
624:     PetscFEGetQuadrature((PetscFE) obj, &quad);
625:     DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
626:     PetscFEGetDimension((PetscFE) obj, &Nb);
627:     PetscFEGetNumComponents((PetscFE) obj, &Nc);
628:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
629:     PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);
630:     PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
631:     DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
632:     DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
633:     PetscArrayzero(gradient, cdim*Nc);
634:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
635:       PetscInt qc = 0, q;

637:       PetscDSGetDiscretization(prob, f, &obj);
638:       PetscArrayzero(interpolant,Nc);
639:       PetscArrayzero(interpolantGrad, cdim*Nc);
640:       for (q = 0; q < Nq; ++q) {
641:         PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, x, &fegeom, q, interpolant, interpolantGrad);
642:         for (fc = 0; fc < Nc; ++fc) {
643:           const PetscReal wt = quadWeights[q*qNc+qc+fc];

645:           field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
646:           for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
647:         }
648:       }
649:       fieldOffset += Nb;
650:       qc          += Nc;
651:     }
652:     PetscFree2(interpolant, interpolantGrad);
653:     DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
654:     for (fc = 0; fc < Nc; ++fc) {
655:       field[fc] /= cg.volume;
656:       for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
657:     }
658:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);
659:     PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
660:   }
661:   return(0);
662: }

664: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
665: {
666:   PetscDS        prob;
667:   void          *ctx;
668:   MPI_Comm       comm;
669:   PetscInt       numAdapt = adaptor->numSeq, adaptIter;
670:   PetscInt       dim, coordDim, numFields, cStart, cEnd, c;

674:   DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
675:   VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
676:   PetscObjectGetComm((PetscObject) adaptor, &comm);
677:   DMGetDimension(adaptor->idm, &dim);
678:   DMGetCoordinateDim(adaptor->idm, &coordDim);
679:   DMGetApplicationContext(adaptor->idm, &ctx);
680:   DMGetDS(adaptor->idm, &prob);
681:   PetscDSGetNumFields(prob, &numFields);

683:   /* Adapt until nothing changes */
684:   /* Adapt for a specified number of iterates */
685:   for (adaptIter = 0; adaptIter < numAdapt-1; ++adaptIter) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));}
686:   for (adaptIter = 0; adaptIter < numAdapt;   ++adaptIter) {
687:     PetscBool adapted = PETSC_FALSE;
688:     DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
689:     Vec       x       = adaptIter ? *ax  : inx, locX, ox;

691:     DMGetLocalVector(dm, &locX);
692:     DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
693:     DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
694:     DMAdaptorPreAdapt(adaptor, locX);
695:     if (doSolve) {
696:       SNES snes;

698:       DMAdaptorGetSolver(adaptor, &snes);
699:       SNESSolve(snes, NULL, adaptIter ? *ax : x);
700:     }
701:     /* DMAdaptorMonitor(adaptor);
702:        Print iterate, memory used, DM, solution */
703:     switch (adaptor->adaptCriterion) {
704:     case DM_ADAPTATION_REFINE:
705:       DMRefine(dm, comm, &odm);
706:       if (!odm) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
707:       adapted = PETSC_TRUE;
708:       break;
709:     case DM_ADAPTATION_LABEL:
710:     {
711:       /* Adapt DM
712:            Create local solution
713:            Reconstruct gradients (FVM) or solve adjoint equation (FEM)
714:            Produce cellwise error indicator */
715:       DM                 plex;
716:       DMLabel            adaptLabel;
717:       IS                 refineIS, coarsenIS;
718:       Vec                errVec;
719:       PetscScalar       *errArray;
720:       const PetscScalar *pointSols;
721:       PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
722:       PetscInt           nRefine, nCoarsen;

724:       DMConvert(dm, DMPLEX, &plex);
725:       DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
726:       DMPlexGetInteriorCellStratum(plex, &cStart, &cEnd);

728:       VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
729:       VecSetUp(errVec);
730:       VecGetArray(errVec, &errArray);
731:       VecGetArrayRead(locX, &pointSols);
732:       for (c = cStart; c < cEnd; ++c) {
733:         PetscReal errInd;

735:         DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);
736:         errArray[c-cStart] = errInd;
737:         minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
738:         minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
739:       }
740:       VecRestoreArrayRead(locX, &pointSols);
741:       VecRestoreArray(errVec, &errArray);
742:       PetscGlobalMinMaxReal(PetscObjectComm((PetscObject) adaptor), minMaxInd, minMaxIndGlobal);
743:       PetscInfo2(adaptor, "DMAdaptor: error indicator range (%E, %E)\n", minMaxIndGlobal[0], minMaxIndGlobal[1]);
744:       /*     Compute IS from VecTagger */
745:       VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS);
746:       VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS);
747:       ISGetSize(refineIS, &nRefine);
748:       ISGetSize(coarsenIS, &nCoarsen);
749:       PetscInfo2(adaptor, "DMAdaptor: numRefine %D, numCoarsen %D\n", nRefine, nCoarsen);
750:       if (nRefine)  {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE,  refineIS);}
751:       if (nCoarsen) {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);}
752:       ISDestroy(&coarsenIS);
753:       ISDestroy(&refineIS);
754:       VecDestroy(&errVec);
755:       /*     Adapt DM from label */
756:       if (nRefine || nCoarsen) {
757:         DMAdaptLabel(dm, adaptLabel, &odm);
758:         adapted = PETSC_TRUE;
759:       }
760:       DMLabelDestroy(&adaptLabel);
761:       DMDestroy(&plex);
762:     }
763:     break;
764:     case DM_ADAPTATION_METRIC:
765:     {
766:       DM           dmGrad,   dmHess,   dmMetric;
767:       PetscDS      probGrad, probHess;
768:       Vec          xGrad,    xHess,    metric;
769:       PetscSection sec, msec;
770:       PetscScalar *H, *M, integral;
771:       PetscReal    N;
772:       DMLabel      bdLabel;
773:       PetscInt     Nd = coordDim*coordDim, f, vStart, vEnd, v;

775:       /*     Compute vertexwise gradients from cellwise gradients */
776:       DMClone(dm, &dmGrad);
777:       DMClone(dm, &dmHess);
778:       DMGetDS(dmGrad, &probGrad);
779:       DMGetDS(dmHess, &probHess);
780:       for (f = 0; f < numFields; ++f) {
781:         PetscFE         fe, feGrad, feHess;
782:         PetscDualSpace  Q;
783:         DM              K;
784:         PetscQuadrature q;
785:         PetscInt        Nc, qorder;
786:         const char     *prefix;

788:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
789:         PetscFEGetNumComponents(fe, &Nc);
790:         PetscFEGetDualSpace(fe, &Q);
791:         PetscDualSpaceGetDM(Q, &K);
792:         DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
793:         PetscFEGetQuadrature(fe, &q);
794:         PetscQuadratureGetOrder(q, &qorder);
795:         PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);
796:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feGrad);
797:         PetscDSSetDiscretization(probGrad, f, (PetscObject) feGrad);
798:         PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feHess);
799:         PetscDSSetDiscretization(probHess, f, (PetscObject) feHess);
800:         PetscFEDestroy(&feGrad);
801:         PetscFEDestroy(&feHess);
802:       }
803:       DMGetGlobalVector(dmGrad, &xGrad);
804:       VecViewFromOptions(x, NULL, "-sol_adapt_loc_pre_view");
805:       DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
806:       VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
807:       /*     Compute vertexwise Hessians from cellwise Hessians */
808:       DMGetGlobalVector(dmHess, &xHess);
809:       DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
810:       VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
811:       /*     Compute metric */
812:       DMClone(dm, &dmMetric);
813:       DMGetLocalSection(dm, &sec);
814:       DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
815:       PetscSectionCreate(PetscObjectComm((PetscObject) dm), &msec);
816:       PetscSectionSetNumFields(msec, 1);
817:       PetscSectionSetFieldComponents(msec, 0, Nd);
818:       PetscSectionSetChart(msec, vStart, vEnd);
819:       for (v = vStart; v < vEnd; ++v) {
820:         PetscSectionSetDof(msec, v, Nd);
821:         PetscSectionSetFieldDof(msec, v, 0, Nd);
822:       }
823:       PetscSectionSetUp(msec);
824:       DMSetLocalSection(dmMetric, msec);
825:       PetscSectionDestroy(&msec);
826:       DMGetLocalVector(dmMetric, &metric);
827:       /*       N is the target size */
828:       N    = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
829:       if (adaptor->monitor) {PetscPrintf(PETSC_COMM_SELF, "N_orig: %D N_adapt: %g\n", vEnd - vStart, N);}
830:       /*       |H| means take the absolute value of eigenvalues */
831:       VecGetArray(xHess, &H);
832:       VecGetArray(metric, &M);
833:       for (v = vStart; v < vEnd; ++v) {
834:         PetscScalar *Hp;

836:         DMPlexPointLocalRef(dmHess, v, H, &Hp);
837:         DMAdaptorModifyHessian_Private(coordDim, adaptor->h_min, adaptor->h_max, Hp);
838:       }
839:       /*       Pointwise on vertices M(x) = N^{2/d} (\int_\Omega det(|H|)^{p/(2p+d)})^{-2/d} det(|H|)^{-1/(2p+d)} |H| for L_p */
840:       PetscDSSetObjective(probHess, 0, detHFunc);
841:       DMPlexComputeIntegralFEM(dmHess, xHess, &integral, NULL);
842:       for (v = vStart; v < vEnd; ++v) {
843:         const PetscInt     p = 1;
844:         const PetscScalar *Hp;
845:         PetscScalar       *Mp;
846:         PetscReal          detH, fact;
847:         PetscInt           i;

849:         DMPlexPointLocalRead(dmHess, v, H, (void *) &Hp);
850:         DMPlexPointLocalRef(dmMetric, v, M, &Mp);
851:         if      (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, Hp);
852:         else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, Hp);
853:         else SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "Dimension %d not supported", dim);
854:         fact = PetscPowReal(N, 2.0/dim) * PetscPowReal(PetscRealPart(integral), -2.0/dim) * PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim));
855: #if 0
856:         PetscPrintf(PETSC_COMM_SELF, "fact: %g integral: %g |detH|: %g termA: %g termB: %g\n", fact, integral, PetscAbsReal(detH), PetscPowReal(integral, -2.0/dim), PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim)));
857:         DMPrintCellMatrix(v, "H", coordDim, coordDim, Hp);
858: #endif
859:         for (i = 0; i < Nd; ++i) {
860:           Mp[i] = fact * Hp[i];
861:         }
862:       }
863:       VecRestoreArray(xHess, &H);
864:       VecRestoreArray(metric, &M);
865:       /*     Adapt DM from metric */
866:       DMGetLabel(dm, "marker", &bdLabel);
867:       DMAdaptMetric(dm, metric, bdLabel, &odm);
868:       adapted = PETSC_TRUE;
869:       /* Cleanup */
870:       DMRestoreLocalVector(dmMetric, &metric);
871:       DMDestroy(&dmMetric);
872:       DMRestoreGlobalVector(dmHess, &xHess);
873:       DMRestoreGlobalVector(dmGrad, &xGrad);
874:       DMDestroy(&dmGrad);
875:       DMDestroy(&dmHess);
876:     }
877:     break;
878:     default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
879:     }
880:     DMAdaptorPostAdapt(adaptor);
881:     DMRestoreLocalVector(dm, &locX);
882:     /* If DM was adapted, replace objects and recreate solution */
883:     if (adapted) {
884:       const char *name;

886:       PetscObjectGetName((PetscObject) dm, &name);
887:       PetscObjectSetName((PetscObject) odm, name);
888:       /* Reconfigure solver */
889:       SNESReset(adaptor->snes);
890:       SNESSetDM(adaptor->snes, odm);
891:       DMAdaptorSetSolver(adaptor, adaptor->snes);
892:       DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
893:       SNESSetFromOptions(adaptor->snes);
894:       /* Transfer system */
895:       DMCopyDisc(adaptor->idm, odm);
896:       /* Transfer solution */
897:       DMCreateGlobalVector(odm, &ox);
898:       PetscObjectGetName((PetscObject) x, &name);
899:       PetscObjectSetName((PetscObject) ox, name);
900:       DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
901:       /* Cleanup adaptivity info */
902:       if (adaptIter > 0) {PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));}
903:       DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
904:       DMDestroy(&dm);
905:       VecDestroy(&x);
906:       *adm = odm;
907:       *ax  = ox;
908:     } else {
909:       *adm = dm;
910:       *ax  = x;
911:       adaptIter = numAdapt;
912:     }
913:     if (adaptIter < numAdapt-1) {
914:       DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
915:       VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
916:     }
917:   }
918:   DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
919:   VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
920:   return(0);
921: }

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

926:   Not collective

928:   Input Parameter:
929: + adaptor  - The DMAdaptor object
930: . x        - The global approximate solution
931: - strategy - The adaptation strategy

933:   Output Parameters:
934: + adm - The adapted DM
935: - ax  - The adapted solution

937:   Options database keys:
938: + -snes_adapt <strategy> : initial, sequential, multigrid
939: . -adapt_gradient_view : View the Clement interpolant of the solution gradient
940: . -adapt_hessian_view : View the Clement interpolant of the solution Hessian
941: - -adapt_metric_view : View the metric tensor for adaptive mesh refinement

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

948:   Level: intermediate

950: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
951: @*/
952: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
953: {

957:   switch (strategy)
958:   {
959:   case DM_ADAPTATION_INITIAL:
960:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
961:     break;
962:   case DM_ADAPTATION_SEQUENTIAL:
963:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
964:     break;
965:   default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
966:   }
967:   return(0);
968: }