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: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`, `PetscConvEst`, `PetscConvEstCreate()`
34: @*/
35: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
36: {
37: VecTaggerBox refineBox, coarsenBox;
39: PetscFunctionBegin;
40: PetscAssertPointer(adaptor, 2);
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: [](ch_dmbase), `DM`, `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 values in 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: [](ch_dmbase), `DM`, `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: [](ch_dmbase), `DM`, `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: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
161: @*/
162: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
163: {
164: PetscFunctionBegin;
166: PetscAssertPointer(snes, 2);
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, this MUST have an attached `DM`/`PetscDS`, so that the exact solution can be computed
180: Level: intermediate
182: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
183: @*/
184: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
185: {
186: PetscFunctionBegin;
189: adaptor->snes = snes;
190: PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*@
195: DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter
197: Not Collective
199: Input Parameter:
200: . adaptor - The `DMAdaptor` object
202: Output Parameter:
203: . num - The number of adaptations
205: Level: intermediate
207: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
208: @*/
209: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
210: {
211: PetscFunctionBegin;
213: PetscAssertPointer(num, 2);
214: *num = adaptor->numSeq;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@
219: DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
221: Not Collective
223: Input Parameters:
224: + adaptor - The `DMAdaptor` object
225: - num - The number of adaptations
227: Level: intermediate
229: .seealso: [](ch_dmbase), `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
230: @*/
231: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
232: {
233: PetscFunctionBegin;
235: adaptor->numSeq = num;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: DMAdaptorSetUp - After the solver is specified, creates data structures for controlling adaptivity
242: Collective
244: Input Parameter:
245: . adaptor - The `DMAdaptor` object
247: Level: beginner
249: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
250: @*/
251: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
252: {
253: PetscDS prob;
254: PetscInt Nf, f;
256: PetscFunctionBegin;
257: PetscCall(DMGetDS(adaptor->idm, &prob));
258: PetscCall(VecTaggerSetUp(adaptor->refineTag));
259: PetscCall(VecTaggerSetUp(adaptor->coarsenTag));
260: PetscCall(PetscDSGetNumFields(prob, &Nf));
261: PetscCall(PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx));
262: for (f = 0; f < Nf; ++f) {
263: PetscCall(PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]));
264: /* TODO Have a flag that forces projection rather than using the exact solution */
265: if (adaptor->exactSol[0]) PetscCall(DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private));
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
271: {
272: PetscFunctionBegin;
273: *tfunc = adaptor->ops->transfersolution;
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
278: {
279: PetscFunctionBegin;
280: adaptor->ops->transfersolution = tfunc;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
285: {
286: DM plex;
287: PetscDS prob;
288: PetscObject obj;
289: PetscClassId id;
290: PetscBool isForest;
292: PetscFunctionBegin;
293: PetscCall(DMConvert(adaptor->idm, DMPLEX, &plex));
294: PetscCall(DMGetDS(adaptor->idm, &prob));
295: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
296: PetscCall(PetscObjectGetClassId(obj, &id));
297: PetscCall(DMIsForest(adaptor->idm, &isForest));
298: if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
299: if (isForest) {
300: adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
301: }
302: #if defined(PETSC_HAVE_PRAGMATIC)
303: else {
304: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
305: }
306: #elif defined(PETSC_HAVE_MMG)
307: else {
308: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
309: }
310: #elif defined(PETSC_HAVE_PARMMG)
311: else {
312: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
313: }
314: #else
315: else {
316: adaptor->adaptCriterion = DM_ADAPTATION_REFINE;
317: }
318: #endif
319: }
320: if (id == PETSCFV_CLASSID) {
321: adaptor->femType = PETSC_FALSE;
322: } else {
323: adaptor->femType = PETSC_TRUE;
324: }
325: if (adaptor->femType) {
326: /* Compute local solution bc */
327: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
328: } else {
329: PetscFV fvm = (PetscFV)obj;
330: PetscLimiter noneLimiter;
331: Vec grad;
333: PetscCall(PetscFVGetComputeGradients(fvm, &adaptor->computeGradient));
334: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
335: /* Use no limiting when reconstructing gradients for adaptivity */
336: PetscCall(PetscFVGetLimiter(fvm, &adaptor->limiter));
337: PetscCall(PetscObjectReference((PetscObject)adaptor->limiter));
338: PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
339: PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
340: PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
341: /* Get FVM data */
342: PetscCall(DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM));
343: PetscCall(VecGetDM(adaptor->cellGeom, &adaptor->cellDM));
344: PetscCall(VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
345: /* Compute local solution bc */
346: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
347: /* Compute gradients */
348: PetscCall(DMCreateGlobalVector(adaptor->gradDM, &grad));
349: PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
350: PetscCall(DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad));
351: PetscCall(DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
352: PetscCall(DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
353: PetscCall(VecDestroy(&grad));
354: PetscCall(VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
355: }
356: PetscCall(DMDestroy(&plex));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
361: {
362: PetscReal time = 0.0;
363: Mat interp;
364: void *ctx;
366: PetscFunctionBegin;
367: PetscCall(DMGetApplicationContext(dm, &ctx));
368: if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
369: else {
370: switch (adaptor->adaptCriterion) {
371: case DM_ADAPTATION_LABEL:
372: PetscCall(DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time));
373: break;
374: case DM_ADAPTATION_REFINE:
375: case DM_ADAPTATION_METRIC:
376: PetscCall(DMCreateInterpolation(dm, adm, &interp, NULL));
377: PetscCall(MatInterpolate(interp, x, ax));
378: PetscCall(DMInterpolate(dm, interp, adm));
379: PetscCall(MatDestroy(&interp));
380: break;
381: default:
382: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
383: }
384: }
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
389: {
390: PetscDS prob;
391: PetscObject obj;
392: PetscClassId id;
394: PetscFunctionBegin;
395: PetscCall(DMGetDS(adaptor->idm, &prob));
396: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
397: PetscCall(PetscObjectGetClassId(obj, &id));
398: if (id == PETSCFV_CLASSID) {
399: PetscFV fvm = (PetscFV)obj;
401: PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
402: /* Restore original limiter */
403: PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));
405: PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
406: PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
407: PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
408: }
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*
413: DMAdaptorSimpleErrorIndicator - Use the integrated gradient as an error indicator in the `DMAdaptor`
415: Input Parameters:
416: + adaptor - The `DMAdaptor` object
417: . dim - The topological dimension
418: . cell - The cell
419: . field - The field integrated over the cell
420: . gradient - The gradient integrated over the cell
421: . cg - A `PetscFVCellGeom` struct
422: - ctx - A user context
424: Output Parameter:
425: . errInd - The error indicator
427: Developer Note:
428: Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API
430: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorComputeErrorIndicator()`
431: */
432: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
433: {
434: PetscReal err = 0.;
435: PetscInt c, d;
437: PetscFunctionBeginHot;
438: for (c = 0; c < Nc; c++) {
439: for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
440: }
441: *errInd = cg->volume * err;
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
446: {
447: PetscDS prob;
448: PetscObject obj;
449: PetscClassId id;
450: void *ctx;
451: PetscQuadrature quad;
452: PetscInt dim, d, cdim, Nc;
454: PetscFunctionBegin;
455: *errInd = 0.;
456: PetscCall(DMGetDimension(plex, &dim));
457: PetscCall(DMGetCoordinateDim(plex, &cdim));
458: PetscCall(DMGetApplicationContext(plex, &ctx));
459: PetscCall(DMGetDS(plex, &prob));
460: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
461: PetscCall(PetscObjectGetClassId(obj, &id));
462: if (id == PETSCFV_CLASSID) {
463: const PetscScalar *pointSols;
464: const PetscScalar *pointSol;
465: const PetscScalar *pointGrad;
466: PetscFVCellGeom *cg;
468: PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
469: PetscCall(VecGetArrayRead(locX, &pointSols));
470: PetscCall(DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol));
471: PetscCall(DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad));
472: PetscCall(DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg));
473: PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
474: PetscCall(VecRestoreArrayRead(locX, &pointSols));
475: } else {
476: PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
477: PetscFVCellGeom cg;
478: PetscFEGeom fegeom;
479: const PetscReal *quadWeights;
480: PetscReal *coords;
481: PetscInt Nb, fc, Nq, qNc, Nf, f, fieldOffset;
483: fegeom.dim = dim;
484: fegeom.dimEmbed = cdim;
485: PetscCall(PetscDSGetNumFields(prob, &Nf));
486: PetscCall(PetscFEGetQuadrature((PetscFE)obj, &quad));
487: PetscCall(DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x));
488: PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
489: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
490: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights));
491: PetscCall(PetscCalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ));
492: PetscCall(PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad));
493: PetscCall(DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
494: PetscCall(DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL));
495: PetscCall(PetscArrayzero(gradient, cdim * Nc));
496: for (f = 0, fieldOffset = 0; f < Nf; ++f) {
497: PetscInt qc = 0, q;
499: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
500: PetscCall(PetscArrayzero(interpolant, Nc));
501: PetscCall(PetscArrayzero(interpolantGrad, cdim * Nc));
502: for (q = 0; q < Nq; ++q) {
503: PetscCall(PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad));
504: for (fc = 0; fc < Nc; ++fc) {
505: const PetscReal wt = quadWeights[q * qNc + qc + fc];
507: field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
508: for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
509: }
510: }
511: fieldOffset += Nb;
512: qc += Nc;
513: }
514: PetscCall(PetscFree2(interpolant, interpolantGrad));
515: PetscCall(DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x));
516: for (fc = 0; fc < Nc; ++fc) {
517: field[fc] /= cg.volume;
518: for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
519: }
520: PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, field, gradient, &cg, errInd, ctx);
521: PetscCall(PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
522: }
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: 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[])
527: {
528: PetscInt i, j;
530: for (i = 0; i < dim; ++i) {
531: for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
532: }
533: }
535: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
536: {
537: PetscDS prob;
538: void *ctx;
539: MPI_Comm comm;
540: PetscInt numAdapt = adaptor->numSeq, adaptIter;
541: PetscInt dim, coordDim, numFields, cStart, cEnd, c;
543: PetscFunctionBegin;
544: PetscCall(DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view"));
545: PetscCall(VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view"));
546: PetscCall(PetscObjectGetComm((PetscObject)adaptor, &comm));
547: PetscCall(DMGetDimension(adaptor->idm, &dim));
548: PetscCall(DMGetCoordinateDim(adaptor->idm, &coordDim));
549: PetscCall(DMGetApplicationContext(adaptor->idm, &ctx));
550: PetscCall(DMGetDS(adaptor->idm, &prob));
551: PetscCall(PetscDSGetNumFields(prob, &numFields));
552: PetscCheck(numFields != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fields is zero!");
554: /* Adapt until nothing changes */
555: /* Adapt for a specified number of iterates */
556: for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm)));
557: for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
558: PetscBool adapted = PETSC_FALSE;
559: DM dm = adaptIter ? *adm : adaptor->idm, odm;
560: Vec x = adaptIter ? *ax : inx, locX, ox;
562: PetscCall(DMGetLocalVector(dm, &locX));
563: PetscCall(DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
564: PetscCall(DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
565: PetscCall(DMAdaptorPreAdapt(adaptor, locX));
566: if (doSolve) {
567: SNES snes;
569: PetscCall(DMAdaptorGetSolver(adaptor, &snes));
570: PetscCall(SNESSolve(snes, NULL, adaptIter ? *ax : x));
571: }
572: /* PetscCall(DMAdaptorMonitor(adaptor));
573: Print iterate, memory used, DM, solution */
574: switch (adaptor->adaptCriterion) {
575: case DM_ADAPTATION_REFINE:
576: PetscCall(DMRefine(dm, comm, &odm));
577: PetscCheck(odm, comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
578: adapted = PETSC_TRUE;
579: break;
580: case DM_ADAPTATION_LABEL: {
581: /* Adapt DM
582: Create local solution
583: Reconstruct gradients (FVM) or solve adjoint equation (FEM)
584: Produce cellwise error indicator */
585: DM plex;
586: DMLabel adaptLabel;
587: IS refineIS, coarsenIS;
588: Vec errVec;
589: PetscScalar *errArray;
590: const PetscScalar *pointSols;
591: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
592: PetscInt nRefine, nCoarsen;
594: PetscCall(DMConvert(dm, DMPLEX, &plex));
595: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
596: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
598: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)adaptor), cEnd - cStart, PETSC_DETERMINE, &errVec));
599: PetscCall(VecSetUp(errVec));
600: PetscCall(VecGetArray(errVec, &errArray));
601: PetscCall(VecGetArrayRead(locX, &pointSols));
602: for (c = cStart; c < cEnd; ++c) {
603: PetscReal errInd;
605: PetscCall(DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd));
606: errArray[c - cStart] = errInd;
607: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
608: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
609: }
610: PetscCall(VecRestoreArrayRead(locX, &pointSols));
611: PetscCall(VecRestoreArray(errVec, &errArray));
612: PetscCall(PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxIndGlobal));
613: PetscCall(PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxIndGlobal[0], (double)minMaxIndGlobal[1]));
614: /* Compute IS from VecTagger */
615: PetscCall(VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS, NULL));
616: PetscCall(VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS, NULL));
617: PetscCall(ISGetSize(refineIS, &nRefine));
618: PetscCall(ISGetSize(coarsenIS, &nCoarsen));
619: PetscCall(PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen));
620: if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
621: if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
622: PetscCall(ISDestroy(&coarsenIS));
623: PetscCall(ISDestroy(&refineIS));
624: PetscCall(VecDestroy(&errVec));
625: /* Adapt DM from label */
626: if (nRefine || nCoarsen) {
627: PetscCall(DMAdaptLabel(dm, adaptLabel, &odm));
628: adapted = PETSC_TRUE;
629: }
630: PetscCall(DMLabelDestroy(&adaptLabel));
631: PetscCall(DMDestroy(&plex));
632: } break;
633: case DM_ADAPTATION_METRIC: {
634: DM dmGrad, dmHess, dmMetric, dmDet;
635: Vec xGrad, xHess, metric, determinant;
636: PetscReal N;
637: DMLabel bdLabel = NULL, rgLabel = NULL;
638: PetscBool higherOrder = PETSC_FALSE;
639: PetscInt Nd = coordDim * coordDim, f, vStart, vEnd;
640: 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[]);
642: PetscCall(PetscMalloc(1, &funcs));
643: funcs[0] = identityFunc;
645: /* Setup finite element spaces */
646: PetscCall(DMClone(dm, &dmGrad));
647: PetscCall(DMClone(dm, &dmHess));
648: PetscCheck(numFields <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple fields not yet considered"); // TODO
649: for (f = 0; f < numFields; ++f) {
650: PetscFE fe, feGrad, feHess;
651: PetscDualSpace Q;
652: PetscSpace space;
653: DM K;
654: PetscQuadrature q;
655: PetscInt Nc, qorder, p;
656: const char *prefix;
658: PetscCall(PetscDSGetDiscretization(prob, f, (PetscObject *)&fe));
659: PetscCall(PetscFEGetNumComponents(fe, &Nc));
660: PetscCheck(Nc <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple components not yet considered"); // TODO
661: PetscCall(PetscFEGetBasisSpace(fe, &space));
662: PetscCall(PetscSpaceGetDegree(space, NULL, &p));
663: if (p > 1) higherOrder = PETSC_TRUE;
664: PetscCall(PetscFEGetDualSpace(fe, &Q));
665: PetscCall(PetscDualSpaceGetDM(Q, &K));
666: PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
667: PetscCall(PetscFEGetQuadrature(fe, &q));
668: PetscCall(PetscQuadratureGetOrder(q, &qorder));
669: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix));
670: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad));
671: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess));
672: PetscCall(DMSetField(dmGrad, f, NULL, (PetscObject)feGrad));
673: PetscCall(DMSetField(dmHess, f, NULL, (PetscObject)feHess));
674: PetscCall(DMCreateDS(dmGrad));
675: PetscCall(DMCreateDS(dmHess));
676: PetscCall(PetscFEDestroy(&feGrad));
677: PetscCall(PetscFEDestroy(&feHess));
678: }
679: /* Compute vertexwise gradients from cellwise gradients */
680: PetscCall(DMCreateLocalVector(dmGrad, &xGrad));
681: PetscCall(VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view"));
682: PetscCall(DMPlexComputeGradientClementInterpolant(dm, locX, xGrad));
683: PetscCall(VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view"));
684: /* Compute vertexwise Hessians from cellwise Hessians */
685: PetscCall(DMCreateLocalVector(dmHess, &xHess));
686: PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess));
687: PetscCall(VecViewFromOptions(xHess, NULL, "-adapt_hessian_view"));
688: PetscCall(VecDestroy(&xGrad));
689: PetscCall(DMDestroy(&dmGrad));
690: /* Compute L-p normalized metric */
691: PetscCall(DMClone(dm, &dmMetric));
692: N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
693: if (adaptor->monitor) {
694: PetscMPIInt rank, size;
695: PetscCallMPI(MPI_Comm_rank(comm, &size));
696: PetscCallMPI(MPI_Comm_rank(comm, &rank));
697: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] N_orig: %" PetscInt_FMT " N_adapt: %g\n", rank, vEnd - vStart, (double)N));
698: }
699: PetscCall(DMPlexMetricSetTargetComplexity(dmMetric, (PetscReal)N));
700: if (higherOrder) {
701: /* Project Hessian into P1 space, if required */
702: PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
703: PetscCall(DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric));
704: PetscCall(VecDestroy(&xHess));
705: xHess = metric;
706: }
707: PetscCall(PetscFree(funcs));
708: PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
709: PetscCall(DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet));
710: PetscCall(DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant));
711: PetscCall(VecDestroy(&determinant));
712: PetscCall(DMDestroy(&dmDet));
713: PetscCall(VecDestroy(&xHess));
714: PetscCall(DMDestroy(&dmHess));
715: /* Adapt DM from metric */
716: PetscCall(DMGetLabel(dm, "marker", &bdLabel));
717: PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm));
718: adapted = PETSC_TRUE;
719: /* Cleanup */
720: PetscCall(VecDestroy(&metric));
721: PetscCall(DMDestroy(&dmMetric));
722: } break;
723: default:
724: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
725: }
726: PetscCall(DMAdaptorPostAdapt(adaptor));
727: PetscCall(DMRestoreLocalVector(dm, &locX));
728: /* If DM was adapted, replace objects and recreate solution */
729: if (adapted) {
730: const char *name;
732: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
733: PetscCall(PetscObjectSetName((PetscObject)odm, name));
734: /* Reconfigure solver */
735: PetscCall(SNESReset(adaptor->snes));
736: PetscCall(SNESSetDM(adaptor->snes, odm));
737: PetscCall(DMAdaptorSetSolver(adaptor, adaptor->snes));
738: PetscCall(DMPlexSetSNESLocalFEM(odm, PETSC_FALSE, ctx));
739: PetscCall(SNESSetFromOptions(adaptor->snes));
740: /* Transfer system */
741: PetscCall(DMCopyDisc(dm, odm));
742: /* Transfer solution */
743: PetscCall(DMCreateGlobalVector(odm, &ox));
744: PetscCall(PetscObjectGetName((PetscObject)x, &name));
745: PetscCall(PetscObjectSetName((PetscObject)ox, name));
746: PetscCall(DMAdaptorTransferSolution(adaptor, dm, x, odm, ox));
747: /* Cleanup adaptivity info */
748: if (adaptIter > 0) PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm)));
749: PetscCall(DMForestSetAdaptivityForest(dm, NULL)); /* clear internal references to the previous dm */
750: PetscCall(DMDestroy(&dm));
751: PetscCall(VecDestroy(&x));
752: *adm = odm;
753: *ax = ox;
754: } else {
755: *adm = dm;
756: *ax = x;
757: adaptIter = numAdapt;
758: }
759: if (adaptIter < numAdapt - 1) {
760: PetscCall(DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view"));
761: PetscCall(VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view"));
762: }
763: }
764: PetscCall(DMViewFromOptions(*adm, NULL, "-dm_adapt_view"));
765: PetscCall(VecViewFromOptions(*ax, NULL, "-sol_adapt_view"));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /*@
770: DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem
772: Not Collective
774: Input Parameters:
775: + adaptor - The `DMAdaptor` object
776: . x - The global approximate solution
777: - strategy - The adaptation strategy, see `DMAdaptationStrategy`
779: Output Parameters:
780: + adm - The adapted `DM`
781: - ax - The adapted solution
783: Options Database Keys:
784: + -snes_adapt <strategy> - initial, sequential, multigrid
785: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
786: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
787: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement
789: Level: intermediate
791: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptationStrategy`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`
792: @*/
793: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
794: {
795: PetscFunctionBegin;
796: switch (strategy) {
797: case DM_ADAPTATION_INITIAL:
798: PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
799: break;
800: case DM_ADAPTATION_SEQUENTIAL:
801: PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
802: break;
803: default:
804: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
805: }
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }