Actual source code: dmadapt.c
1: #include <petscdmadaptor.h>
2: #include <petscdmplex.h>
3: #include <petscdmforest.h>
4: #include <petscds.h>
5: #include <petscblaslapack.h>
7: #include <petsc/private/dmadaptorimpl.h>
8: #include <petsc/private/dmpleximpl.h>
9: #include <petsc/private/petscfeimpl.h>
11: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor, PetscInt, PetscInt, const PetscScalar *, const PetscScalar *, const PetscFVCellGeom *, PetscReal *, void *);
13: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
14: {
16: DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);
17: return 0;
18: }
20: /*@
21: DMAdaptorCreate - Create a DMAdaptor object. Its purpose is to construct a adaptation DMLabel or metric Vec that can be used to modify the DM.
23: Collective
25: Input Parameter:
26: . comm - The communicator for the DMAdaptor object
28: Output Parameter:
29: . adaptor - The DMAdaptor object
31: Level: beginner
33: .seealso: DMAdaptorDestroy(), DMAdaptorAdapt()
34: @*/
35: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
36: {
37: VecTaggerBox refineBox, coarsenBox;
40: PetscSysInitializePackage();
41: PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);
43: (*adaptor)->monitor = PETSC_FALSE;
44: (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
45: (*adaptor)->numSeq = 1;
46: (*adaptor)->Nadapt = -1;
47: (*adaptor)->refinementFactor = 2.0;
48: (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
49: refineBox.min = refineBox.max = PETSC_MAX_REAL;
50: VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);
51: PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");
52: VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);
53: VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);
54: coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
55: VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);
56: PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");
57: VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);
58: VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);
59: return 0;
60: }
62: /*@
63: DMAdaptorDestroy - Destroys a DMAdaptor object
65: Collective on DMAdaptor
67: Input Parameter:
68: . adaptor - The DMAdaptor object
70: Level: beginner
72: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
73: @*/
74: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
75: {
76: if (!*adaptor) return 0;
78: if (--((PetscObject)(*adaptor))->refct > 0) {
79: *adaptor = NULL;
80: return 0;
81: }
82: VecTaggerDestroy(&(*adaptor)->refineTag);
83: VecTaggerDestroy(&(*adaptor)->coarsenTag);
84: PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);
85: PetscHeaderDestroy(adaptor);
86: return 0;
87: }
89: /*@
90: DMAdaptorSetFromOptions - Sets a DMAdaptor object from options
92: Collective on DMAdaptor
94: Input Parameter:
95: . adaptor - The DMAdaptor object
97: Options Database Keys:
98: + -adaptor_monitor <bool> - Monitor the adaptation process
99: . -adaptor_sequence_num <num> - Number of adaptations to generate an optimal grid
100: . -adaptor_target_num <num> - Set the target number of vertices N_adapt, -1 for automatic determination
101: - -adaptor_refinement_factor <r> - Set r such that N_adapt = r^dim N_orig
103: Level: beginner
105: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
106: @*/
107: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
108: {
111: PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");
112: PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
113: PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
114: PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
115: PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
116: PetscOptionsEnd();
117: VecTaggerSetFromOptions(adaptor->refineTag);
118: VecTaggerSetFromOptions(adaptor->coarsenTag);
119: return 0;
120: }
122: /*@
123: DMAdaptorView - Views a DMAdaptor object
125: Collective on DMAdaptor
127: Input Parameters:
128: + adaptor - The DMAdaptor object
129: - viewer - The PetscViewer object
131: Level: beginner
133: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
134: @*/
135: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
136: {
137: PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
138: PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
139: PetscViewerASCIIPrintf(viewer, " sequence length: %D\n", adaptor->numSeq);
140: VecTaggerView(adaptor->refineTag, viewer);
141: VecTaggerView(adaptor->coarsenTag, viewer);
142: return 0;
143: }
145: /*@
146: DMAdaptorGetSolver - Gets the solver used to produce discrete solutions
148: Not collective
150: Input Parameter:
151: . adaptor - The DMAdaptor object
153: Output Parameter:
154: . snes - The solver
156: Level: intermediate
158: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
159: @*/
160: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
161: {
164: *snes = adaptor->snes;
165: return 0;
166: }
168: /*@
169: DMAdaptorSetSolver - Sets the solver used to produce discrete solutions
171: Not collective
173: Input Parameters:
174: + adaptor - The DMAdaptor object
175: - snes - The solver
177: Level: intermediate
179: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
181: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
182: @*/
183: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
184: {
187: adaptor->snes = snes;
188: SNESGetDM(adaptor->snes, &adaptor->idm);
189: return 0;
190: }
192: /*@
193: DMAdaptorGetSequenceLength - Gets the number of sequential adaptations
195: Not collective
197: Input Parameter:
198: . adaptor - The DMAdaptor object
200: Output Parameter:
201: . num - The number of adaptations
203: Level: intermediate
205: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
206: @*/
207: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
208: {
211: *num = adaptor->numSeq;
212: return 0;
213: }
215: /*@
216: DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
218: Not collective
220: Input Parameters:
221: + adaptor - The DMAdaptor object
222: - num - The number of adaptations
224: Level: intermediate
226: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
227: @*/
228: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
229: {
231: adaptor->numSeq = num;
232: return 0;
233: }
235: /*@
236: DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity
238: Collective on DMAdaptor
240: Input Parameters:
241: . adaptor - The DMAdaptor object
243: Level: beginner
245: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
246: @*/
247: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
248: {
249: PetscDS prob;
250: PetscInt Nf, f;
252: DMGetDS(adaptor->idm, &prob);
253: VecTaggerSetUp(adaptor->refineTag);
254: VecTaggerSetUp(adaptor->coarsenTag);
255: PetscDSGetNumFields(prob, &Nf);
256: PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
257: for (f = 0; f < Nf; ++f) {
258: PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
259: /* TODO Have a flag that forces projection rather than using the exact solution */
260: if (adaptor->exactSol[0]) DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);
261: }
262: return 0;
263: }
265: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
266: {
267: *tfunc = adaptor->ops->transfersolution;
268: return 0;
269: }
271: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
272: {
273: adaptor->ops->transfersolution = tfunc;
274: return 0;
275: }
277: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
278: {
279: DM plex;
280: PetscDS prob;
281: PetscObject obj;
282: PetscClassId id;
283: PetscBool isForest;
285: DMConvert(adaptor->idm, DMPLEX, &plex);
286: DMGetDS(adaptor->idm, &prob);
287: PetscDSGetDiscretization(prob, 0, &obj);
288: PetscObjectGetClassId(obj, &id);
289: DMIsForest(adaptor->idm, &isForest);
290: if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
291: if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
292: #if defined(PETSC_HAVE_PRAGMATIC)
293: else {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
294: #elif defined(PETSC_HAVE_MMG)
295: else {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
296: #elif defined(PETSC_HAVE_PARMMG)
297: else {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
298: #else
299: else {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
300: #endif
301: }
302: if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
303: else {adaptor->femType = PETSC_TRUE;}
304: if (adaptor->femType) {
305: /* Compute local solution bc */
306: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
307: } else {
308: PetscFV fvm = (PetscFV) obj;
309: PetscLimiter noneLimiter;
310: Vec grad;
312: PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
313: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
314: /* Use no limiting when reconstructing gradients for adaptivity */
315: PetscFVGetLimiter(fvm, &adaptor->limiter);
316: PetscObjectReference((PetscObject) adaptor->limiter);
317: PetscLimiterCreate(PetscObjectComm((PetscObject) fvm), &noneLimiter);
318: PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
319: PetscFVSetLimiter(fvm, noneLimiter);
320: /* Get FVM data */
321: DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
322: VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
323: VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
324: /* Compute local solution bc */
325: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
326: /* Compute gradients */
327: DMCreateGlobalVector(adaptor->gradDM, &grad);
328: DMPlexReconstructGradientsFVM(plex, locX, grad);
329: DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
330: DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
331: DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
332: VecDestroy(&grad);
333: VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
334: }
335: DMDestroy(&plex);
336: return 0;
337: }
339: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
340: {
341: PetscReal time = 0.0;
342: Mat interp;
343: void *ctx;
345: DMGetApplicationContext(dm, &ctx);
346: if (adaptor->ops->transfersolution) {
347: (*adaptor->ops->transfersolution)(adaptor, dm, x, adm, ax, ctx);
348: } else {
349: switch (adaptor->adaptCriterion) {
350: case DM_ADAPTATION_LABEL:
351: DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
352: break;
353: case DM_ADAPTATION_REFINE:
354: case DM_ADAPTATION_METRIC:
355: DMCreateInterpolation(dm, adm, &interp, NULL);
356: MatInterpolate(interp, x, ax);
357: DMInterpolate(dm, interp, adm);
358: MatDestroy(&interp);
359: break;
360: default: SETERRQ(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %D", adaptor->adaptCriterion);
361: }
362: }
363: return 0;
364: }
366: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
367: {
368: PetscDS prob;
369: PetscObject obj;
370: PetscClassId id;
372: DMGetDS(adaptor->idm, &prob);
373: PetscDSGetDiscretization(prob, 0, &obj);
374: PetscObjectGetClassId(obj, &id);
375: if (id == PETSCFV_CLASSID) {
376: PetscFV fvm = (PetscFV) obj;
378: PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
379: /* Restore original limiter */
380: PetscFVSetLimiter(fvm, adaptor->limiter);
382: VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
383: VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
384: DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
385: }
386: return 0;
387: }
389: /*
390: DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator
392: Input Parameters:
393: + adaptor - The DMAdaptor object
394: . dim - The topological dimension
395: . cell - The cell
396: . field - The field integrated over the cell
397: . gradient - The gradient integrated over the cell
398: . cg - A PetscFVCellGeom struct
399: - ctx - A user context
401: Output Parameter:
402: . errInd - The error indicator
404: .seealso: DMAdaptorComputeErrorIndicator()
405: */
406: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
407: {
408: PetscReal err = 0.;
409: PetscInt c, d;
412: for (c = 0; c < Nc; c++) {
413: for (d = 0; d < dim; ++d) {
414: err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
415: }
416: }
417: *errInd = cg->volume * err;
418: return 0;
419: }
421: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
422: {
423: PetscDS prob;
424: PetscObject obj;
425: PetscClassId id;
426: void *ctx;
427: PetscQuadrature quad;
428: PetscInt dim, d, cdim, Nc;
430: *errInd = 0.;
431: DMGetDimension(plex, &dim);
432: DMGetCoordinateDim(plex, &cdim);
433: DMGetApplicationContext(plex, &ctx);
434: DMGetDS(plex, &prob);
435: PetscDSGetDiscretization(prob, 0, &obj);
436: PetscObjectGetClassId(obj, &id);
437: if (id == PETSCFV_CLASSID) {
438: const PetscScalar *pointSols;
439: const PetscScalar *pointSol;
440: const PetscScalar *pointGrad;
441: PetscFVCellGeom *cg;
443: PetscFVGetNumComponents((PetscFV) obj, &Nc);
444: VecGetArrayRead(locX, &pointSols);
445: DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
446: DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
447: DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
448: (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
449: VecRestoreArrayRead(locX, &pointSols);
450: } else {
451: PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
452: PetscFVCellGeom cg;
453: PetscFEGeom fegeom;
454: const PetscReal *quadWeights;
455: PetscReal *coords;
456: PetscInt Nb, fc, Nq, qNc, Nf, f, fieldOffset;
458: fegeom.dim = dim;
459: fegeom.dimEmbed = cdim;
460: PetscDSGetNumFields(prob, &Nf);
461: PetscFEGetQuadrature((PetscFE) obj, &quad);
462: DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
463: PetscFEGetDimension((PetscFE) obj, &Nb);
464: PetscFEGetNumComponents((PetscFE) obj, &Nc);
465: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
466: PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);
467: PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
468: DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
469: DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
470: PetscArrayzero(gradient, cdim*Nc);
471: for (f = 0, fieldOffset = 0; f < Nf; ++f) {
472: PetscInt qc = 0, q;
474: PetscDSGetDiscretization(prob, f, &obj);
475: PetscArrayzero(interpolant,Nc);
476: PetscArrayzero(interpolantGrad, cdim*Nc);
477: for (q = 0; q < Nq; ++q) {
478: PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, 1, x, &fegeom, q, interpolant, interpolantGrad);
479: for (fc = 0; fc < Nc; ++fc) {
480: const PetscReal wt = quadWeights[q*qNc+qc+fc];
482: field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
483: for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
484: }
485: }
486: fieldOffset += Nb;
487: qc += Nc;
488: }
489: PetscFree2(interpolant, interpolantGrad);
490: DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
491: for (fc = 0; fc < Nc; ++fc) {
492: field[fc] /= cg.volume;
493: for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
494: }
495: (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);
496: PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
497: }
498: return 0;
499: }
501: static void identityFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
502: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
503: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
504: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
505: {
506: PetscInt i, j;
508: for (i = 0; i < dim; ++i) {
509: for (j = 0; j < dim; ++j) {
510: f[i+dim*j] = u[i+dim*j];
511: }
512: }
513: }
515: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
516: {
517: PetscDS prob;
518: void *ctx;
519: MPI_Comm comm;
520: PetscInt numAdapt = adaptor->numSeq, adaptIter;
521: PetscInt dim, coordDim, numFields, cStart, cEnd, c;
523: DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
524: VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
525: PetscObjectGetComm((PetscObject) adaptor, &comm);
526: DMGetDimension(adaptor->idm, &dim);
527: DMGetCoordinateDim(adaptor->idm, &coordDim);
528: DMGetApplicationContext(adaptor->idm, &ctx);
529: DMGetDS(adaptor->idm, &prob);
530: PetscDSGetNumFields(prob, &numFields);
533: /* Adapt until nothing changes */
534: /* Adapt for a specified number of iterates */
535: for (adaptIter = 0; adaptIter < numAdapt-1; ++adaptIter) PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));
536: for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
537: PetscBool adapted = PETSC_FALSE;
538: DM dm = adaptIter ? *adm : adaptor->idm, odm;
539: Vec x = adaptIter ? *ax : inx, locX, ox;
541: DMGetLocalVector(dm, &locX);
542: DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
543: DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
544: DMAdaptorPreAdapt(adaptor, locX);
545: if (doSolve) {
546: SNES snes;
548: DMAdaptorGetSolver(adaptor, &snes);
549: SNESSolve(snes, NULL, adaptIter ? *ax : x);
550: }
551: /* DMAdaptorMonitor(adaptor);
552: Print iterate, memory used, DM, solution */
553: switch (adaptor->adaptCriterion) {
554: case DM_ADAPTATION_REFINE:
555: DMRefine(dm, comm, &odm);
557: adapted = PETSC_TRUE;
558: break;
559: case DM_ADAPTATION_LABEL:
560: {
561: /* Adapt DM
562: Create local solution
563: Reconstruct gradients (FVM) or solve adjoint equation (FEM)
564: Produce cellwise error indicator */
565: DM plex;
566: DMLabel adaptLabel;
567: IS refineIS, coarsenIS;
568: Vec errVec;
569: PetscScalar *errArray;
570: const PetscScalar *pointSols;
571: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
572: PetscInt nRefine, nCoarsen;
574: DMConvert(dm, DMPLEX, &plex);
575: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
576: DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);
578: VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
579: VecSetUp(errVec);
580: VecGetArray(errVec, &errArray);
581: VecGetArrayRead(locX, &pointSols);
582: for (c = cStart; c < cEnd; ++c) {
583: PetscReal errInd;
585: DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);
586: errArray[c-cStart] = errInd;
587: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
588: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
589: }
590: VecRestoreArrayRead(locX, &pointSols);
591: VecRestoreArray(errVec, &errArray);
592: PetscGlobalMinMaxReal(PetscObjectComm((PetscObject) adaptor), minMaxInd, minMaxIndGlobal);
593: PetscInfo(adaptor, "DMAdaptor: error indicator range (%E, %E)\n", minMaxIndGlobal[0], minMaxIndGlobal[1]);
594: /* Compute IS from VecTagger */
595: VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS,NULL);
596: VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS,NULL);
597: ISGetSize(refineIS, &nRefine);
598: ISGetSize(coarsenIS, &nCoarsen);
599: PetscInfo(adaptor, "DMAdaptor: numRefine %D, numCoarsen %D\n", nRefine, nCoarsen);
600: if (nRefine) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS);
601: if (nCoarsen) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);
602: ISDestroy(&coarsenIS);
603: ISDestroy(&refineIS);
604: VecDestroy(&errVec);
605: /* Adapt DM from label */
606: if (nRefine || nCoarsen) {
607: DMAdaptLabel(dm, adaptLabel, &odm);
608: adapted = PETSC_TRUE;
609: }
610: DMLabelDestroy(&adaptLabel);
611: DMDestroy(&plex);
612: }
613: break;
614: case DM_ADAPTATION_METRIC:
615: {
616: DM dmGrad, dmHess, dmMetric;
617: Vec xGrad, xHess, metric;
618: PetscReal N;
619: DMLabel bdLabel = NULL, rgLabel = NULL;
620: PetscBool higherOrder = PETSC_FALSE;
621: PetscInt Nd = coordDim*coordDim, f, vStart, vEnd;
622: void (**funcs)(PetscInt, PetscInt, PetscInt,
623: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
624: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
625: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
627: PetscMalloc(1, &funcs);
628: funcs[0] = identityFunc;
630: /* Setup finite element spaces */
631: DMClone(dm, &dmGrad);
632: DMClone(dm, &dmHess);
634: for (f = 0; f < numFields; ++f) {
635: PetscFE fe, feGrad, feHess;
636: PetscDualSpace Q;
637: PetscSpace space;
638: DM K;
639: PetscQuadrature q;
640: PetscInt Nc, qorder, p;
641: const char *prefix;
643: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
644: PetscFEGetNumComponents(fe, &Nc);
646: PetscFEGetBasisSpace(fe, &space);
647: PetscSpaceGetDegree(space, NULL, &p);
648: if (p > 1) higherOrder = PETSC_TRUE;
649: PetscFEGetDualSpace(fe, &Q);
650: PetscDualSpaceGetDM(Q, &K);
651: DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
652: PetscFEGetQuadrature(fe, &q);
653: PetscQuadratureGetOrder(q, &qorder);
654: PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);
655: PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, PETSC_TRUE, prefix, qorder, &feGrad);
656: PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, PETSC_TRUE, prefix, qorder, &feHess);
657: DMSetField(dmGrad, f, NULL, (PetscObject)feGrad);
658: DMSetField(dmHess, f, NULL, (PetscObject)feHess);
659: DMCreateDS(dmGrad);
660: DMCreateDS(dmHess);
661: PetscFEDestroy(&feGrad);
662: PetscFEDestroy(&feHess);
663: }
664: /* Compute vertexwise gradients from cellwise gradients */
665: DMCreateLocalVector(dmGrad, &xGrad);
666: VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view");
667: DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
668: VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
669: /* Compute vertexwise Hessians from cellwise Hessians */
670: DMCreateLocalVector(dmHess, &xHess);
671: DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
672: VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
673: VecDestroy(&xGrad);
674: DMDestroy(&dmGrad);
675: /* Compute L-p normalized metric */
676: DMClone(dm, &dmMetric);
677: N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
678: if (adaptor->monitor) {
679: PetscMPIInt rank, size;
680: MPI_Comm_rank(comm, &size);
681: MPI_Comm_rank(comm, &rank);
682: PetscPrintf(PETSC_COMM_SELF, "[%D] N_orig: %D N_adapt: %g\n", rank, vEnd - vStart, N);
683: }
684: DMPlexMetricSetTargetComplexity(dmMetric, (PetscReal) N);
685: if (higherOrder) {
686: /* Project Hessian into P1 space, if required */
687: DMPlexMetricCreate(dmMetric, 0, &metric);
688: DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric);
689: VecDestroy(&xHess);
690: xHess = metric;
691: }
692: PetscFree(funcs);
693: DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, &metric);
694: VecDestroy(&xHess);
695: DMDestroy(&dmHess);
696: /* Adapt DM from metric */
697: DMGetLabel(dm, "marker", &bdLabel);
698: DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm);
699: adapted = PETSC_TRUE;
700: /* Cleanup */
701: VecDestroy(&metric);
702: DMDestroy(&dmMetric);
703: }
704: break;
705: default: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
706: }
707: DMAdaptorPostAdapt(adaptor);
708: DMRestoreLocalVector(dm, &locX);
709: /* If DM was adapted, replace objects and recreate solution */
710: if (adapted) {
711: const char *name;
713: PetscObjectGetName((PetscObject) dm, &name);
714: PetscObjectSetName((PetscObject) odm, name);
715: /* Reconfigure solver */
716: SNESReset(adaptor->snes);
717: SNESSetDM(adaptor->snes, odm);
718: DMAdaptorSetSolver(adaptor, adaptor->snes);
719: DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
720: SNESSetFromOptions(adaptor->snes);
721: /* Transfer system */
722: DMCopyDisc(dm, odm);
723: /* Transfer solution */
724: DMCreateGlobalVector(odm, &ox);
725: PetscObjectGetName((PetscObject) x, &name);
726: PetscObjectSetName((PetscObject) ox, name);
727: DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
728: /* Cleanup adaptivity info */
729: if (adaptIter > 0) PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));
730: DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
731: DMDestroy(&dm);
732: VecDestroy(&x);
733: *adm = odm;
734: *ax = ox;
735: } else {
736: *adm = dm;
737: *ax = x;
738: adaptIter = numAdapt;
739: }
740: if (adaptIter < numAdapt-1) {
741: DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
742: VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
743: }
744: }
745: DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
746: VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
747: return 0;
748: }
750: /*@
751: DMAdaptorAdapt - Creates a new DM that is adapted to the problem
753: Not collective
755: Input Parameters:
756: + adaptor - The DMAdaptor object
757: . x - The global approximate solution
758: - strategy - The adaptation strategy
760: Output Parameters:
761: + adm - The adapted DM
762: - ax - The adapted solution
764: Options database keys:
765: + -snes_adapt <strategy> - initial, sequential, multigrid
766: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
767: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
768: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement
770: Note: The available adaptation strategies are:
771: $ 1) Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
772: $ 2) Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
773: $ 3) Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid
775: Level: intermediate
777: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
778: @*/
779: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
780: {
781: switch (strategy)
782: {
783: case DM_ADAPTATION_INITIAL:
784: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
785: break;
786: case DM_ADAPTATION_SEQUENTIAL:
787: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
788: break;
789: default: SETERRQ(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
790: }
791: return 0;
792: }