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