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: `DMAdaptor`, `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
67: Input Parameter:
68: . adaptor - The `DMAdaptor` object
70: Level: beginner
72: .seealso: `DMAdaptor`, `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 properties of a `DMAdaptor` object from the options database
92: Collective
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: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
106: @*/
107: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
108: {
109: PetscOptionsBegin(PetscObjectComm((PetscObject)adaptor), "", "DM Adaptor Options", "DMAdaptor");
110: PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
111: PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
112: PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
113: PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
114: PetscOptionsEnd();
115: VecTaggerSetFromOptions(adaptor->refineTag);
116: VecTaggerSetFromOptions(adaptor->coarsenTag);
117: return 0;
118: }
120: /*@
121: DMAdaptorView - Views a `DMAdaptor` object
123: Collective
125: Input Parameters:
126: + adaptor - The `DMAdaptor` object
127: - viewer - The `PetscViewer` object
129: Level: beginner
131: .seealso: `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
132: @*/
133: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
134: {
135: PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer);
136: PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
137: PetscViewerASCIIPrintf(viewer, " sequence length: %" PetscInt_FMT "\n", adaptor->numSeq);
138: VecTaggerView(adaptor->refineTag, viewer);
139: VecTaggerView(adaptor->coarsenTag, viewer);
140: return 0;
141: }
143: /*@
144: DMAdaptorGetSolver - Gets the solver used to produce discrete solutions
146: Not collective
148: Input Parameter:
149: . adaptor - The `DMAdaptor` object
151: Output Parameter:
152: . snes - The solver
154: Level: intermediate
156: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
157: @*/
158: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
159: {
162: *snes = adaptor->snes;
163: return 0;
164: }
166: /*@
167: DMAdaptorSetSolver - Sets the solver used to produce discrete solutions
169: Not collective
171: Input Parameters:
172: + adaptor - The `DMAdaptor` object
173: - snes - The solver
175: Level: intermediate
177: Note:
178: The solver MUST have an attached `DM`/`DS`, so that we know the exact solution
180: .seealso: `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
181: @*/
182: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
183: {
186: adaptor->snes = snes;
187: SNESGetDM(adaptor->snes, &adaptor->idm);
188: return 0;
189: }
191: /*@
192: DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter
194: Not collective
196: Input Parameter:
197: . adaptor - The `DMAdaptor` object
199: Output Parameter:
200: . num - The number of adaptations
202: Level: intermediate
204: .seealso: `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
205: @*/
206: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
207: {
210: *num = adaptor->numSeq;
211: return 0;
212: }
214: /*@
215: DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
217: Not collective
219: Input Parameters:
220: + adaptor - The `DMAdaptor` object
221: - num - The number of adaptations
223: Level: intermediate
225: .seealso: `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
226: @*/
227: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
228: {
230: adaptor->numSeq = num;
231: return 0;
232: }
234: /*@
235: DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity
237: Collective
239: Input Parameters:
240: . adaptor - The `DMAdaptor` object
242: Level: beginner
244: .seealso: `DMAdaptorCreate()`, `DMAdaptorAdapt()`
245: @*/
246: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
247: {
248: PetscDS prob;
249: PetscInt Nf, f;
251: DMGetDS(adaptor->idm, &prob);
252: VecTaggerSetUp(adaptor->refineTag);
253: VecTaggerSetUp(adaptor->coarsenTag);
254: PetscDSGetNumFields(prob, &Nf);
255: PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
256: for (f = 0; f < Nf; ++f) {
257: PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
258: /* TODO Have a flag that forces projection rather than using the exact solution */
259: if (adaptor->exactSol[0]) DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);
260: }
261: return 0;
262: }
264: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
265: {
266: *tfunc = adaptor->ops->transfersolution;
267: return 0;
268: }
270: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
271: {
272: adaptor->ops->transfersolution = tfunc;
273: return 0;
274: }
276: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
277: {
278: DM plex;
279: PetscDS prob;
280: PetscObject obj;
281: PetscClassId id;
282: PetscBool isForest;
284: DMConvert(adaptor->idm, DMPLEX, &plex);
285: DMGetDS(adaptor->idm, &prob);
286: PetscDSGetDiscretization(prob, 0, &obj);
287: PetscObjectGetClassId(obj, &id);
288: DMIsForest(adaptor->idm, &isForest);
289: if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
290: if (isForest) {
291: adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
292: }
293: #if defined(PETSC_HAVE_PRAGMATIC)
294: else {
295: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
296: }
297: #elif defined(PETSC_HAVE_MMG)
298: else {
299: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
300: }
301: #elif defined(PETSC_HAVE_PARMMG)
302: else {
303: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
304: }
305: #else
306: else {
307: adaptor->adaptCriterion = DM_ADAPTATION_REFINE;
308: }
309: #endif
310: }
311: if (id == PETSCFV_CLASSID) {
312: adaptor->femType = PETSC_FALSE;
313: } else {
314: adaptor->femType = PETSC_TRUE;
315: }
316: if (adaptor->femType) {
317: /* Compute local solution bc */
318: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
319: } else {
320: PetscFV fvm = (PetscFV)obj;
321: PetscLimiter noneLimiter;
322: Vec grad;
324: PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
325: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
326: /* Use no limiting when reconstructing gradients for adaptivity */
327: PetscFVGetLimiter(fvm, &adaptor->limiter);
328: PetscObjectReference((PetscObject)adaptor->limiter);
329: PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter);
330: PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
331: PetscFVSetLimiter(fvm, noneLimiter);
332: /* Get FVM data */
333: DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
334: VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
335: VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
336: /* Compute local solution bc */
337: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
338: /* Compute gradients */
339: DMCreateGlobalVector(adaptor->gradDM, &grad);
340: DMPlexReconstructGradientsFVM(plex, locX, grad);
341: DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
342: DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
343: DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
344: VecDestroy(&grad);
345: VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
346: }
347: DMDestroy(&plex);
348: return 0;
349: }
351: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
352: {
353: PetscReal time = 0.0;
354: Mat interp;
355: void *ctx;
357: DMGetApplicationContext(dm, &ctx);
358: if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
359: else {
360: switch (adaptor->adaptCriterion) {
361: case DM_ADAPTATION_LABEL:
362: DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
363: break;
364: case DM_ADAPTATION_REFINE:
365: case DM_ADAPTATION_METRIC:
366: DMCreateInterpolation(dm, adm, &interp, NULL);
367: MatInterpolate(interp, x, ax);
368: DMInterpolate(dm, interp, adm);
369: MatDestroy(&interp);
370: break;
371: default:
372: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
373: }
374: }
375: return 0;
376: }
378: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
379: {
380: PetscDS prob;
381: PetscObject obj;
382: PetscClassId id;
384: DMGetDS(adaptor->idm, &prob);
385: PetscDSGetDiscretization(prob, 0, &obj);
386: PetscObjectGetClassId(obj, &id);
387: if (id == PETSCFV_CLASSID) {
388: PetscFV fvm = (PetscFV)obj;
390: PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
391: /* Restore original limiter */
392: PetscFVSetLimiter(fvm, adaptor->limiter);
394: VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
395: VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
396: DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
397: }
398: return 0;
399: }
401: /*
402: DMAdaptorSimpleErrorIndicator - Use the integrated gradient as an error indicator
404: Input Parameters:
405: + adaptor - The `DMAdaptor` object
406: . dim - The topological dimension
407: . cell - The cell
408: . field - The field integrated over the cell
409: . gradient - The gradient integrated over the cell
410: . cg - A `PetscFVCellGeom` struct
411: - ctx - A user context
413: Output Parameter:
414: . errInd - The error indicator
416: Developer Note:
417: Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API
419: .seealso: `DMAdaptor`, `DMAdaptorComputeErrorIndicator()`
420: */
421: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
422: {
423: PetscReal err = 0.;
424: PetscInt c, d;
427: for (c = 0; c < Nc; c++) {
428: for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
429: }
430: *errInd = cg->volume * err;
431: return 0;
432: }
434: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
435: {
436: PetscDS prob;
437: PetscObject obj;
438: PetscClassId id;
439: void *ctx;
440: PetscQuadrature quad;
441: PetscInt dim, d, cdim, Nc;
443: *errInd = 0.;
444: DMGetDimension(plex, &dim);
445: DMGetCoordinateDim(plex, &cdim);
446: DMGetApplicationContext(plex, &ctx);
447: DMGetDS(plex, &prob);
448: PetscDSGetDiscretization(prob, 0, &obj);
449: PetscObjectGetClassId(obj, &id);
450: if (id == PETSCFV_CLASSID) {
451: const PetscScalar *pointSols;
452: const PetscScalar *pointSol;
453: const PetscScalar *pointGrad;
454: PetscFVCellGeom *cg;
456: PetscFVGetNumComponents((PetscFV)obj, &Nc);
457: VecGetArrayRead(locX, &pointSols);
458: DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol);
459: DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad);
460: DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
461: PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
462: VecRestoreArrayRead(locX, &pointSols);
463: } else {
464: PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
465: PetscFVCellGeom cg;
466: PetscFEGeom fegeom;
467: const PetscReal *quadWeights;
468: PetscReal *coords;
469: PetscInt Nb, fc, Nq, qNc, Nf, f, fieldOffset;
471: fegeom.dim = dim;
472: fegeom.dimEmbed = cdim;
473: PetscDSGetNumFields(prob, &Nf);
474: PetscFEGetQuadrature((PetscFE)obj, &quad);
475: DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
476: PetscFEGetDimension((PetscFE)obj, &Nb);
477: PetscFEGetNumComponents((PetscFE)obj, &Nc);
478: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
479: PetscMalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ);
480: PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad);
481: DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
482: DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
483: PetscArrayzero(gradient, cdim * Nc);
484: for (f = 0, fieldOffset = 0; f < Nf; ++f) {
485: PetscInt qc = 0, q;
487: PetscDSGetDiscretization(prob, f, &obj);
488: PetscArrayzero(interpolant, Nc);
489: PetscArrayzero(interpolantGrad, cdim * Nc);
490: for (q = 0; q < Nq; ++q) {
491: PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad);
492: for (fc = 0; fc < Nc; ++fc) {
493: const PetscReal wt = quadWeights[q * qNc + qc + fc];
495: field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
496: for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
497: }
498: }
499: fieldOffset += Nb;
500: qc += Nc;
501: }
502: PetscFree2(interpolant, interpolantGrad);
503: DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
504: for (fc = 0; fc < Nc; ++fc) {
505: field[fc] /= cg.volume;
506: for (d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
507: }
508: PetscUseTypeMethod(adaptor, computeerrorindicator, dim, Nc, field, gradient, &cg, errInd, ctx);
509: PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ);
510: }
511: return 0;
512: }
514: 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[])
515: {
516: PetscInt i, j;
518: for (i = 0; i < dim; ++i) {
519: for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
520: }
521: }
523: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
524: {
525: PetscDS prob;
526: void *ctx;
527: MPI_Comm comm;
528: PetscInt numAdapt = adaptor->numSeq, adaptIter;
529: PetscInt dim, coordDim, numFields, cStart, cEnd, c;
531: DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
532: VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
533: PetscObjectGetComm((PetscObject)adaptor, &comm);
534: DMGetDimension(adaptor->idm, &dim);
535: DMGetCoordinateDim(adaptor->idm, &coordDim);
536: DMGetApplicationContext(adaptor->idm, &ctx);
537: DMGetDS(adaptor->idm, &prob);
538: PetscDSGetNumFields(prob, &numFields);
541: /* Adapt until nothing changes */
542: /* Adapt for a specified number of iterates */
543: for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));
544: for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
545: PetscBool adapted = PETSC_FALSE;
546: DM dm = adaptIter ? *adm : adaptor->idm, odm;
547: Vec x = adaptIter ? *ax : inx, locX, ox;
549: DMGetLocalVector(dm, &locX);
550: DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
551: DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
552: DMAdaptorPreAdapt(adaptor, locX);
553: if (doSolve) {
554: SNES snes;
556: DMAdaptorGetSolver(adaptor, &snes);
557: SNESSolve(snes, NULL, adaptIter ? *ax : x);
558: }
559: /* DMAdaptorMonitor(adaptor);
560: Print iterate, memory used, DM, solution */
561: switch (adaptor->adaptCriterion) {
562: case DM_ADAPTATION_REFINE:
563: DMRefine(dm, comm, &odm);
565: adapted = PETSC_TRUE;
566: break;
567: case DM_ADAPTATION_LABEL: {
568: /* Adapt DM
569: Create local solution
570: Reconstruct gradients (FVM) or solve adjoint equation (FEM)
571: Produce cellwise error indicator */
572: DM plex;
573: DMLabel adaptLabel;
574: IS refineIS, coarsenIS;
575: Vec errVec;
576: PetscScalar *errArray;
577: const PetscScalar *pointSols;
578: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
579: PetscInt nRefine, nCoarsen;
581: DMConvert(dm, DMPLEX, &plex);
582: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
583: DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);
585: VecCreateMPI(PetscObjectComm((PetscObject)adaptor), cEnd - cStart, PETSC_DETERMINE, &errVec);
586: VecSetUp(errVec);
587: VecGetArray(errVec, &errArray);
588: VecGetArrayRead(locX, &pointSols);
589: for (c = cStart; c < cEnd; ++c) {
590: PetscReal errInd;
592: DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);
593: errArray[c - cStart] = errInd;
594: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
595: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
596: }
597: VecRestoreArrayRead(locX, &pointSols);
598: VecRestoreArray(errVec, &errArray);
599: PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxIndGlobal);
600: PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxIndGlobal[0], (double)minMaxIndGlobal[1]);
601: /* Compute IS from VecTagger */
602: VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS, NULL);
603: VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS, NULL);
604: ISGetSize(refineIS, &nRefine);
605: ISGetSize(coarsenIS, &nCoarsen);
606: PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen);
607: if (nRefine) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS);
608: if (nCoarsen) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);
609: ISDestroy(&coarsenIS);
610: ISDestroy(&refineIS);
611: VecDestroy(&errVec);
612: /* Adapt DM from label */
613: if (nRefine || nCoarsen) {
614: DMAdaptLabel(dm, adaptLabel, &odm);
615: adapted = PETSC_TRUE;
616: }
617: DMLabelDestroy(&adaptLabel);
618: DMDestroy(&plex);
619: } break;
620: case DM_ADAPTATION_METRIC: {
621: DM dmGrad, dmHess, dmMetric, dmDet;
622: Vec xGrad, xHess, metric, determinant;
623: PetscReal N;
624: DMLabel bdLabel = NULL, rgLabel = NULL;
625: PetscBool higherOrder = PETSC_FALSE;
626: PetscInt Nd = coordDim * coordDim, f, vStart, vEnd;
627: 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[]);
629: PetscMalloc(1, &funcs);
630: funcs[0] = identityFunc;
632: /* Setup finite element spaces */
633: DMClone(dm, &dmGrad);
634: DMClone(dm, &dmHess);
636: for (f = 0; f < numFields; ++f) {
637: PetscFE fe, feGrad, feHess;
638: PetscDualSpace Q;
639: PetscSpace space;
640: DM K;
641: PetscQuadrature q;
642: PetscInt Nc, qorder, p;
643: const char *prefix;
645: PetscDSGetDiscretization(prob, f, (PetscObject *)&fe);
646: PetscFEGetNumComponents(fe, &Nc);
648: PetscFEGetBasisSpace(fe, &space);
649: PetscSpaceGetDegree(space, NULL, &p);
650: if (p > 1) higherOrder = PETSC_TRUE;
651: PetscFEGetDualSpace(fe, &Q);
652: PetscDualSpaceGetDM(Q, &K);
653: DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
654: PetscFEGetQuadrature(fe, &q);
655: PetscQuadratureGetOrder(q, &qorder);
656: PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix);
657: PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad);
658: PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess);
659: DMSetField(dmGrad, f, NULL, (PetscObject)feGrad);
660: DMSetField(dmHess, f, NULL, (PetscObject)feHess);
661: DMCreateDS(dmGrad);
662: DMCreateDS(dmHess);
663: PetscFEDestroy(&feGrad);
664: PetscFEDestroy(&feHess);
665: }
666: /* Compute vertexwise gradients from cellwise gradients */
667: DMCreateLocalVector(dmGrad, &xGrad);
668: VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view");
669: DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
670: VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
671: /* Compute vertexwise Hessians from cellwise Hessians */
672: DMCreateLocalVector(dmHess, &xHess);
673: DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
674: VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
675: VecDestroy(&xGrad);
676: DMDestroy(&dmGrad);
677: /* Compute L-p normalized metric */
678: DMClone(dm, &dmMetric);
679: N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
680: if (adaptor->monitor) {
681: PetscMPIInt rank, size;
682: MPI_Comm_rank(comm, &size);
683: MPI_Comm_rank(comm, &rank);
684: PetscPrintf(PETSC_COMM_SELF, "[%d] N_orig: %" PetscInt_FMT " N_adapt: %g\n", rank, vEnd - vStart, (double)N);
685: }
686: DMPlexMetricSetTargetComplexity(dmMetric, (PetscReal)N);
687: if (higherOrder) {
688: /* Project Hessian into P1 space, if required */
689: DMPlexMetricCreate(dmMetric, 0, &metric);
690: DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric);
691: VecDestroy(&xHess);
692: xHess = metric;
693: }
694: PetscFree(funcs);
695: DMPlexMetricCreate(dmMetric, 0, &metric);
696: DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet);
697: DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant);
698: VecDestroy(&determinant);
699: DMDestroy(&dmDet);
700: VecDestroy(&xHess);
701: DMDestroy(&dmHess);
702: /* Adapt DM from metric */
703: DMGetLabel(dm, "marker", &bdLabel);
704: DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm);
705: adapted = PETSC_TRUE;
706: /* Cleanup */
707: VecDestroy(&metric);
708: DMDestroy(&dmMetric);
709: } break;
710: default:
711: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
712: }
713: DMAdaptorPostAdapt(adaptor);
714: DMRestoreLocalVector(dm, &locX);
715: /* If DM was adapted, replace objects and recreate solution */
716: if (adapted) {
717: const char *name;
719: PetscObjectGetName((PetscObject)dm, &name);
720: PetscObjectSetName((PetscObject)odm, name);
721: /* Reconfigure solver */
722: SNESReset(adaptor->snes);
723: SNESSetDM(adaptor->snes, odm);
724: DMAdaptorSetSolver(adaptor, adaptor->snes);
725: DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
726: SNESSetFromOptions(adaptor->snes);
727: /* Transfer system */
728: DMCopyDisc(dm, odm);
729: /* Transfer solution */
730: DMCreateGlobalVector(odm, &ox);
731: PetscObjectGetName((PetscObject)x, &name);
732: PetscObjectSetName((PetscObject)ox, name);
733: DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
734: /* Cleanup adaptivity info */
735: if (adaptIter > 0) PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));
736: DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
737: DMDestroy(&dm);
738: VecDestroy(&x);
739: *adm = odm;
740: *ax = ox;
741: } else {
742: *adm = dm;
743: *ax = x;
744: adaptIter = numAdapt;
745: }
746: if (adaptIter < numAdapt - 1) {
747: DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
748: VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
749: }
750: }
751: DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
752: VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
753: return 0;
754: }
756: /*@
757: DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem
759: Not collective
761: Input Parameters:
762: + adaptor - The `DMAdaptor` object
763: . x - The global approximate solution
764: - strategy - The adaptation strategy
766: Output Parameters:
767: + adm - The adapted `DM`
768: - ax - The adapted solution
770: Options database Keys:
771: + -snes_adapt <strategy> - initial, sequential, multigrid
772: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
773: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
774: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement
776: Note: The available adaptation strategies are:
777: + * - Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
778: . * - Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
779: - * - Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid
781: Level: intermediate
783: .seealso: `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
784: @*/
785: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
786: {
787: switch (strategy) {
788: case DM_ADAPTATION_INITIAL:
789: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
790: break;
791: case DM_ADAPTATION_SEQUENTIAL:
792: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
793: break;
794: default:
795: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
796: }
797: return 0;
798: }