Actual source code: dmadapt.c

petsc-3.11.3 2019-06-26
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, (void **) ctx, 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 on MPI_Comm

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

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

 33:   Level: beginner

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

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

 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)->h_min = 1.;
 54:   (*adaptor)->h_max = 10000.;
 55:   (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
 56:   refineBox.min = refineBox.max = PETSC_MAX_REAL;
 57:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);
 58:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");
 59:   VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);
 60:   VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);
 61:   coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
 62:   VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);
 63:   PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");
 64:   VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);
 65:   VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);
 66:   return(0);
 67: }

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

 72:   Collective on DMAdaptor

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

 77:   Level: beginner

 79: .keywords: DMAdaptor, convergence, destroy
 80: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
 81: @*/
 82: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
 83: {

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

100: /*@
101:   DMAdaptorSetFromOptions - Sets a DMAdaptor object from options

103:   Collective on DMAdaptor

105:   Input Parameters:
106: . adaptor - The DMAdaptor object

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

116:   Level: beginner

118: .keywords: DMAdaptor, convergence, options
119: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
120: @*/
121: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
122: {

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

139: /*@
140:   DMAdaptorView - Views a DMAdaptor object

142:   Collective on DMAdaptor

144:   Input Parameters:
145: + adaptor     - The DMAdaptor object
146: - viewer - The PetscViewer object

148:   Level: beginner

150: .keywords: DMAdaptor, adaptivity, view
151: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
152: @*/
153: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
154: {

158:   PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
159:   PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
160:   PetscViewerASCIIPrintf(viewer, "  sequence length: %D\n", adaptor->numSeq);
161:   VecTaggerView(adaptor->refineTag,  viewer);
162:   VecTaggerView(adaptor->coarsenTag, viewer);
163:   return(0);
164: }

166: /*@
167:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

169:   Not collective

171:   Input Parameter:
172: . adaptor   - The DMAdaptor object

174:   Output Parameter:
175: . snes - The solver

177:   Level: intermediate

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

191: /*@
192:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

194:   Not collective

196:   Input Parameters:
197: + adaptor   - The DMAdaptor object
198: - snes - The solver

200:   Level: intermediate

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

204: .keywords: DMAdaptor, convergence
205: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
206: @*/
207: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
208: {

214:   adaptor->snes = snes;
215:   SNESGetDM(adaptor->snes, &adaptor->idm);
216:   return(0);
217: }

219: /*@
220:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations

222:   Not collective

224:   Input Parameter:
225: . adaptor - The DMAdaptor object

227:   Output Parameter:
228: . num - The number of adaptations

230:   Level: intermediate

232: .keywords: DMAdaptor, convergence
233: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
234: @*/
235: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
236: {
240:   *num = adaptor->numSeq;
241:   return(0);
242: }

244: /*@
245:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

247:   Not collective

249:   Input Parameters:
250: + adaptor - The DMAdaptor object
251: - num - The number of adaptations

253:   Level: intermediate

255: .keywords: DMAdaptor, convergence
256: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
257: @*/
258: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
259: {
262:   adaptor->numSeq = num;
263:   return(0);
264: }

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

269:   Collective on DMAdaptor

271:   Input Parameters:
272: . adaptor - The DMAdaptor object

274:   Level: beginner

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

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

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

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

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

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

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

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

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

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

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

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

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

427: static PetscErrorCode DMAdaptorModifyHessian_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscScalar Hp[])
428: {
429:   PetscScalar   *Hpos;
430:   PetscReal     *eigs;
431:   PetscInt       i, j, k;

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

473:     lwork = 5*dim;
474:     PetscMalloc1(5*dim, &work);
475: #if defined(PETSC_MISSING_LAPACK_GEEV)
476:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
477: #else
478:     {
479:       PetscBLASInt lierr;
480:       PetscBLASInt nb;

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

546: static void detHFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
547:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
548:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
549:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
550: {
551:   const PetscInt p = 1;
552:   PetscReal      detH = 0.0;

554:   if      (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, u);
555:   else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, u);
556:   f0[0] = PetscPowReal(detH, p/(2.*p + dim));
557: }

559: /*
560:   DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator

562:   Input Parameters:
563: + adaptor  - The DMAdaptor object
564: . dim      - The topological dimension
565: . cell     - The cell
566: . field    - The field integrated over the cell
567: . gradient - The gradient integrated over the cell
568: . cg       - A PetscFVCellGeom struct
569: - ctx      - A user context

571:   Output Parameter:
572: . errInd   - The error indicator

574: .seealso: DMAdaptorComputeErrorIndicator()
575: */
576: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
577: {
578:   PetscReal err = 0.;
579:   PetscInt  c, d;

582:   for (c = 0; c < Nc; c++) {
583:     for (d = 0; d < dim; ++d) {
584:       err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
585:     }
586:   }
587:   *errInd = cg->volume * err;
588:   return(0);
589: }

591: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
592: {
593:   PetscDS         prob;
594:   PetscObject     obj;
595:   PetscClassId    id;
596:   void           *ctx;
597:   PetscQuadrature quad;
598:   PetscInt        dim, d, cdim, Nc;
599:   PetscErrorCode  ierr;

602:   *errInd = 0.;
603:   DMGetDimension(plex, &dim);
604:   DMGetCoordinateDim(plex, &cdim);
605:   DMGetApplicationContext(plex, &ctx);
606:   DMGetDS(plex, &prob);
607:   PetscDSGetDiscretization(prob, 0, &obj);
608:   PetscObjectGetClassId(obj, &id);
609:   if (id == PETSCFV_CLASSID) {
610:     const PetscScalar *pointSols;
611:     const PetscScalar *pointSol;
612:     const PetscScalar *pointGrad;
613:     PetscFVCellGeom   *cg;

615:     PetscFVGetNumComponents((PetscFV) obj, &Nc);
616:     VecGetArrayRead(locX, &pointSols);
617:     DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
618:     DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
619:     DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
620:     (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
621:     VecRestoreArrayRead(locX, &pointSols);
622:   } else {
623:     PetscScalar     *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
624:     PetscFVCellGeom  cg;
625:     const PetscReal *quadWeights;
626:     PetscReal       *coords, *detJ, *J, *invJ;
627:     PetscInt         Nb, fc, Nq, qNc, Nf, f, fieldOffset;

629:     PetscDSGetNumFields(prob, &Nf);
630:     PetscFEGetQuadrature((PetscFE) obj, &quad);
631:     DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
632:     PetscFEGetDimension((PetscFE) obj, &Nb);
633:     PetscFEGetNumComponents((PetscFE) obj, &Nc);
634:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
635:     PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&detJ,cdim*cdim*Nq,&J,cdim*cdim*Nq,&invJ);
636:     PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
637:     DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, J, invJ, detJ);
638:     DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
639:     PetscMemzero(gradient, cdim*Nc * sizeof(PetscScalar));
640:     for (f = 0, fieldOffset = 0; f < Nf; ++f) {
641:       PetscInt qc = 0, q;

643:       PetscDSGetDiscretization(prob, f, &obj);
644:       PetscMemzero(interpolant,     Nc * sizeof(PetscScalar));
645:       PetscMemzero(interpolantGrad, cdim*Nc * sizeof(PetscScalar));
646:       for (q = 0; q < Nq; ++q) {
647:         PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, x, dim, invJ, q, interpolant, interpolantGrad);
648:         for (fc = 0; fc < Nc; ++fc) {
649:           const PetscReal wt = quadWeights[q*qNc+qc+fc];

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

670: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
671: {
672:   PetscDS        prob;
673:   void          *ctx;
674:   MPI_Comm       comm;
675:   PetscInt       numAdapt = adaptor->numSeq, adaptIter;
676:   PetscInt       dim, coordDim, numFields, cStart, cEnd, cEndInterior, c;

680:   DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
681:   VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
682:   PetscObjectGetComm((PetscObject) adaptor, &comm);
683:   DMGetDimension(adaptor->idm, &dim);
684:   DMGetCoordinateDim(adaptor->idm, &coordDim);
685:   DMGetApplicationContext(adaptor->idm, &ctx);
686:   DMGetDS(adaptor->idm, &prob);
687:   PetscDSGetNumFields(prob, &numFields);

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

697:     DMGetLocalVector(dm, &locX);
698:     DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
699:     DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
700:     DMAdaptorPreAdapt(adaptor, locX);
701:     if (doSolve) {
702:       SNES snes;

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

730:       DMConvert(dm, DMPLEX, &plex);
731:       DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
732:       DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
733:       DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
734:       cEnd = (cEndInterior < 0) ? cEnd : cEndInterior;

736:       VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
737:       VecSetUp(errVec);
738:       VecGetArray(errVec, &errArray);
739:       VecGetArrayRead(locX, &pointSols);
740:       for (c = cStart; c < cEnd; ++c) {
741:         PetscReal errInd;

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

783:       /*     Compute vertexwise gradients from cellwise gradients */
784:       DMClone(dm, &dmGrad);
785:       DMClone(dm, &dmHess);
786:       DMGetDS(dmGrad, &probGrad);
787:       DMGetDS(dmHess, &probHess);
788:       for (f = 0; f < numFields; ++f) {
789:         PetscFE         fe, feGrad, feHess;
790:         PetscDualSpace  Q;
791:         DM              K;
792:         PetscQuadrature q;
793:         PetscInt        Nc, qorder;
794:         const char     *prefix;

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

844:         DMPlexPointLocalRef(dmHess, v, H, &Hp);
845:         DMAdaptorModifyHessian_Private(coordDim, adaptor->h_min, adaptor->h_max, Hp);
846:       }
847:       /*       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 */
848:       PetscDSSetObjective(probHess, 0, detHFunc);
849:       DMPlexComputeIntegralFEM(dmHess, xHess, &integral, NULL);
850:       for (v = vStart; v < vEnd; ++v) {
851:         const PetscInt     p = 1;
852:         const PetscScalar *Hp;
853:         PetscScalar       *Mp;
854:         PetscReal          detH, fact;
855:         PetscInt           i;

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

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

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

934:   Not collective

936:   Input Parameter:
937: + adaptor  - The DMAdaptor object
938: . x        - The global approximate solution
939: - strategy - The adaptation strategy

941:   Output Parameters:
942: + adm - The adapted DM
943: - ax  - The adapted solution

945:   Options database keys:
946: . -snes_adapt <strategy> : initial, sequential, multigrid
947: . -adapt_gradient_view : View the Clement interpolant of the solution gradient
948: . -adapt_hessian_view : View the Clement interpolant of the solution Hessian
949: . -adapt_metric_view : View the metric tensor for adaptive mesh refinement

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

956:   Level: intermediate

958: .keywords: DMAdaptor, convergence
959: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
960: @*/
961: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
962: {

966:   switch (strategy)
967:   {
968:   case DM_ADAPTATION_INITIAL:
969:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
970:     break;
971:   case DM_ADAPTATION_SEQUENTIAL:
972:     DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
973:     break;
974:   default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
975:   }
976:   return(0);
977: }