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