Actual source code: dmadapt.c
petsc-3.12.5 2020-03-29
1: #include <petscdmadaptor.h>
2: #include <petscdmplex.h>
3: #include <petscdmforest.h>
4: #include <petscds.h>
5: #include <petscblaslapack.h>
7: #include <petsc/private/dmadaptorimpl.h>
8: #include <petsc/private/dmpleximpl.h>
9: #include <petsc/private/petscfeimpl.h>
11: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor, PetscInt, PetscInt, const PetscScalar *, const PetscScalar *, const PetscFVCellGeom *, PetscReal *, void *);
13: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
14: {
18: DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);
19: return(0);
20: }
22: /*@
23: DMAdaptorCreate - Create a DMAdaptor object. Its purpose is to construct a adaptation DMLabel or metric Vec that can be used to modify the DM.
25: Collective
27: Input Parameter:
28: . comm - The communicator for the DMAdaptor object
30: Output Parameter:
31: . adaptor - The DMAdaptor object
33: Level: beginner
35: .seealso: DMAdaptorDestroy(), DMAdaptorAdapt()
36: @*/
37: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
38: {
39: VecTaggerBox refineBox, coarsenBox;
44: PetscSysInitializePackage();
45: PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);
47: (*adaptor)->monitor = PETSC_FALSE;
48: (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
49: (*adaptor)->numSeq = 1;
50: (*adaptor)->Nadapt = -1;
51: (*adaptor)->refinementFactor = 2.0;
52: (*adaptor)->h_min = 1.;
53: (*adaptor)->h_max = 10000.;
54: (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
55: refineBox.min = refineBox.max = PETSC_MAX_REAL;
56: VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);
57: PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");
58: VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);
59: VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);
60: coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
61: VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);
62: PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");
63: VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);
64: VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);
65: return(0);
66: }
68: /*@
69: DMAdaptorDestroy - Destroys a DMAdaptor object
71: Collective on DMAdaptor
73: Input Parameter:
74: . adaptor - The DMAdaptor object
76: Level: beginner
78: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
79: @*/
80: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
81: {
85: if (!*adaptor) return(0);
87: if (--((PetscObject)(*adaptor))->refct > 0) {
88: *adaptor = NULL;
89: return(0);
90: }
91: VecTaggerDestroy(&(*adaptor)->refineTag);
92: VecTaggerDestroy(&(*adaptor)->coarsenTag);
93: PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);
94: PetscHeaderDestroy(adaptor);
95: return(0);
96: }
98: /*@
99: DMAdaptorSetFromOptions - Sets a DMAdaptor object from options
101: Collective on DMAdaptor
103: Input Parameters:
104: . adaptor - The DMAdaptor object
106: Options Database Keys:
107: + -adaptor_monitor <bool> : Monitor the adaptation process
108: . -adaptor_sequence_num <num> : Number of adaptations to generate an optimal grid
109: . -adaptor_target_num <num> : Set the target number of vertices N_adapt, -1 for automatic determination
110: . -adaptor_refinement_factor <r> : Set r such that N_adapt = r^dim N_orig
111: . -adaptor_metric_h_min <min> : Set the minimum eigenvalue of Hessian (sqr max edge length)
112: - -adaptor_metric_h_max <max> : Set the maximum eigenvalue of Hessian (sqr min edge length)
114: Level: beginner
116: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
117: @*/
118: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
119: {
123: PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");
124: PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);
125: PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);
126: PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);
127: PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);
128: PetscOptionsReal("-adaptor_metric_h_min", "Set the minimum eigenvalue of Hessian (sqr max edge length)", "DMAdaptor", adaptor->h_min, &adaptor->h_min, NULL);
129: PetscOptionsReal("-adaptor_metric_h_max", "Set the maximum eigenvalue of Hessian (sqr min edge length)", "DMAdaptor", adaptor->h_max, &adaptor->h_max, NULL);
130: PetscOptionsEnd();
131: VecTaggerSetFromOptions(adaptor->refineTag);
132: VecTaggerSetFromOptions(adaptor->coarsenTag);
133: return(0);
134: }
136: /*@
137: DMAdaptorView - Views a DMAdaptor object
139: Collective on DMAdaptor
141: Input Parameters:
142: + adaptor - The DMAdaptor object
143: - viewer - The PetscViewer object
145: Level: beginner
147: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
148: @*/
149: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
150: {
154: PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);
155: PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");
156: PetscViewerASCIIPrintf(viewer, " sequence length: %D\n", adaptor->numSeq);
157: VecTaggerView(adaptor->refineTag, viewer);
158: VecTaggerView(adaptor->coarsenTag, viewer);
159: return(0);
160: }
162: /*@
163: DMAdaptorGetSolver - Gets the solver used to produce discrete solutions
165: Not collective
167: Input Parameter:
168: . adaptor - The DMAdaptor object
170: Output Parameter:
171: . snes - The solver
173: Level: intermediate
175: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
176: @*/
177: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
178: {
182: *snes = adaptor->snes;
183: return(0);
184: }
186: /*@
187: DMAdaptorSetSolver - Sets the solver used to produce discrete solutions
189: Not collective
191: Input Parameters:
192: + adaptor - The DMAdaptor object
193: - snes - The solver
195: Level: intermediate
197: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
199: .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
200: @*/
201: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
202: {
208: adaptor->snes = snes;
209: SNESGetDM(adaptor->snes, &adaptor->idm);
210: return(0);
211: }
213: /*@
214: DMAdaptorGetSequenceLength - Gets the number of sequential adaptations
216: Not collective
218: Input Parameter:
219: . adaptor - The DMAdaptor object
221: Output Parameter:
222: . num - The number of adaptations
224: Level: intermediate
226: .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
227: @*/
228: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
229: {
233: *num = adaptor->numSeq;
234: return(0);
235: }
237: /*@
238: DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
240: Not collective
242: Input Parameters:
243: + adaptor - The DMAdaptor object
244: - num - The number of adaptations
246: Level: intermediate
248: .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
249: @*/
250: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
251: {
254: adaptor->numSeq = num;
255: return(0);
256: }
258: /*@
259: DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity
261: Collective on DMAdaptor
263: Input Parameters:
264: . adaptor - The DMAdaptor object
266: Level: beginner
268: .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
269: @*/
270: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
271: {
272: PetscDS prob;
273: PetscInt Nf, f;
277: DMGetDS(adaptor->idm, &prob);
278: VecTaggerSetUp(adaptor->refineTag);
279: VecTaggerSetUp(adaptor->coarsenTag);
280: PetscDSGetNumFields(prob, &Nf);
281: PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);
282: for (f = 0; f < Nf; ++f) {
283: PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);
284: /* TODO Have a flag that forces projection rather than using the exact solution */
285: if (adaptor->exactSol[0]) {DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);}
286: }
287: return(0);
288: }
290: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
291: {
293: *tfunc = adaptor->ops->transfersolution;
294: return(0);
295: }
297: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
298: {
300: adaptor->ops->transfersolution = tfunc;
301: return(0);
302: }
304: PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
305: {
306: DM plex;
307: PetscDS prob;
308: PetscObject obj;
309: PetscClassId id;
310: PetscBool isForest;
314: DMConvert(adaptor->idm, DMPLEX, &plex);
315: DMGetDS(adaptor->idm, &prob);
316: PetscDSGetDiscretization(prob, 0, &obj);
317: PetscObjectGetClassId(obj, &id);
318: DMIsForest(adaptor->idm, &isForest);
319: if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
320: if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
321: #if defined(PETSC_HAVE_PRAGMATIC)
322: else {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
323: #else
324: else {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
325: #endif
326: }
327: if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
328: else {adaptor->femType = PETSC_TRUE;}
329: if (adaptor->femType) {
330: /* Compute local solution bc */
331: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
332: } else {
333: PetscFV fvm = (PetscFV) obj;
334: PetscLimiter noneLimiter;
335: Vec grad;
337: PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);
338: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
339: /* Use no limiting when reconstructing gradients for adaptivity */
340: PetscFVGetLimiter(fvm, &adaptor->limiter);
341: PetscObjectReference((PetscObject) adaptor->limiter);
342: PetscLimiterCreate(PetscObjectComm((PetscObject) fvm), &noneLimiter);
343: PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);
344: PetscFVSetLimiter(fvm, noneLimiter);
345: /* Get FVM data */
346: DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);
347: VecGetDM(adaptor->cellGeom, &adaptor->cellDM);
348: VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
349: /* Compute local solution bc */
350: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);
351: /* Compute gradients */
352: DMCreateGlobalVector(adaptor->gradDM, &grad);
353: DMPlexReconstructGradientsFVM(plex, locX, grad);
354: DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);
355: DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
356: DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);
357: VecDestroy(&grad);
358: VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
359: }
360: DMDestroy(&plex);
361: return(0);
362: }
364: PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
365: {
366: PetscReal time = 0.0;
367: Mat interp;
368: void *ctx;
372: DMGetApplicationContext(dm, &ctx);
373: if (adaptor->ops->transfersolution) {
374: (*adaptor->ops->transfersolution)(adaptor, dm, x, adm, ax, ctx);
375: } else {
376: switch (adaptor->adaptCriterion) {
377: case DM_ADAPTATION_LABEL:
378: DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);
379: break;
380: case DM_ADAPTATION_REFINE:
381: case DM_ADAPTATION_METRIC:
382: DMCreateInterpolation(dm, adm, &interp, NULL);
383: MatInterpolate(interp, x, ax);
384: DMInterpolate(dm, interp, adm);
385: MatDestroy(&interp);
386: break;
387: default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %D", adaptor->adaptCriterion);
388: }
389: }
390: return(0);
391: }
393: PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
394: {
395: PetscDS prob;
396: PetscObject obj;
397: PetscClassId id;
401: DMGetDS(adaptor->idm, &prob);
402: PetscDSGetDiscretization(prob, 0, &obj);
403: PetscObjectGetClassId(obj, &id);
404: if (id == PETSCFV_CLASSID) {
405: PetscFV fvm = (PetscFV) obj;
407: PetscFVSetComputeGradients(fvm, adaptor->computeGradient);
408: /* Restore original limiter */
409: PetscFVSetLimiter(fvm, adaptor->limiter);
411: VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);
412: VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);
413: DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);
414: }
415: return(0);
416: }
418: static PetscErrorCode DMAdaptorModifyHessian_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscScalar Hp[])
419: {
420: PetscScalar *Hpos;
421: PetscReal *eigs;
422: PetscInt i, j, k;
426: PetscMalloc2(dim*dim, &Hpos, dim, &eigs);
427: #if 0
428: PetscPrintf(PETSC_COMM_SELF, "H = [");
429: for (i = 0; i < dim; ++i) {
430: if (i > 0) {PetscPrintf(PETSC_COMM_SELF, " ");}
431: for (j = 0; j < dim; ++j) {
432: if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
433: PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
434: }
435: PetscPrintf(PETSC_COMM_SELF, "\n");
436: }
437: PetscPrintf(PETSC_COMM_SELF, "]\n");
438: #endif
439: /* Symmetrize */
440: for (i = 0; i < dim; ++i) {
441: Hpos[i*dim+i] = Hp[i*dim+i];
442: for (j = i+1; j < dim; ++j) {
443: Hpos[i*dim+j] = 0.5*(Hp[i*dim+j] + Hp[j*dim+i]);
444: Hpos[j*dim+i] = Hpos[i*dim+j];
445: }
446: }
447: #if 0
448: PetscPrintf(PETSC_COMM_SELF, "Hs = [");
449: for (i = 0; i < dim; ++i) {
450: if (i > 0) {PetscPrintf(PETSC_COMM_SELF, " ");}
451: for (j = 0; j < dim; ++j) {
452: if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
453: PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
454: }
455: PetscPrintf(PETSC_COMM_SELF, "\n");
456: }
457: PetscPrintf(PETSC_COMM_SELF, "]\n");
458: #endif
459: /* Compute eigendecomposition */
460: {
461: PetscScalar *work;
462: PetscBLASInt lwork;
464: lwork = 5*dim;
465: PetscMalloc1(5*dim, &work);
466: #if defined(PETSC_MISSING_LAPACK_GEEV)
467: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
468: #else
469: {
470: PetscBLASInt lierr;
471: PetscBLASInt nb;
473: PetscBLASIntCast(dim, &nb);
474: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
475: #if defined(PETSC_USE_COMPLEX)
476: {
477: PetscReal *rwork;
478: PetscMalloc1(3*dim, &rwork);
479: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,rwork,&lierr));
480: PetscFree(rwork);
481: }
482: #else
483: PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,&lierr));
484: #endif
485: if (lierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
486: PetscFPTrapPop();
487: }
488: #endif
489: PetscFree(work);
490: }
491: #if 0
492: PetscPrintf(PETSC_COMM_SELF, "L = [");
493: for (i = 0; i < dim; ++i) {
494: if (i > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
495: PetscPrintf(PETSC_COMM_SELF, "%g", eigs[i]);
496: }
497: PetscPrintf(PETSC_COMM_SELF, "]\n");
498: PetscPrintf(PETSC_COMM_SELF, "Q = [");
499: for (i = 0; i < dim; ++i) {
500: if (i > 0) {PetscPrintf(PETSC_COMM_SELF, " ");}
501: for (j = 0; j < dim; ++j) {
502: if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
503: PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);
504: }
505: PetscPrintf(PETSC_COMM_SELF, "\n");
506: }
507: PetscPrintf(PETSC_COMM_SELF, "]\n");
508: #endif
509: /* Reflect to positive orthant, enforce maximum and minimum size, \lambda \propto 1/h^2
510: TODO get domain bounding box */
511: for (i = 0; i < dim; ++i) eigs[i] = PetscMin(h_max, PetscMax(h_min, PetscAbsReal(eigs[i])));
512: /* Reconstruct Hessian */
513: for (i = 0; i < dim; ++i) {
514: for (j = 0; j < dim; ++j) {
515: Hp[i*dim+j] = 0.0;
516: for (k = 0; k < dim; ++k) {
517: Hp[i*dim+j] += Hpos[k*dim+i] * eigs[k] * Hpos[k*dim+j];
518: }
519: }
520: }
521: PetscFree2(Hpos, eigs);
522: #if 0
523: PetscPrintf(PETSC_COMM_SELF, "H+ = [");
524: for (i = 0; i < dim; ++i) {
525: if (i > 0) {PetscPrintf(PETSC_COMM_SELF, " ");}
526: for (j = 0; j < dim; ++j) {
527: if (j > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
528: PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);
529: }
530: PetscPrintf(PETSC_COMM_SELF, "\n");
531: }
532: PetscPrintf(PETSC_COMM_SELF, "]\n");
533: #endif
534: return(0);
535: }
537: static void detHFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
538: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
539: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
540: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
541: {
542: const PetscInt p = 1;
543: PetscReal detH = 0.0;
545: if (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, u);
546: else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, u);
547: f0[0] = PetscPowReal(detH, p/(2.*p + dim));
548: }
550: /*
551: DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator
553: Input Parameters:
554: + adaptor - The DMAdaptor object
555: . dim - The topological dimension
556: . cell - The cell
557: . field - The field integrated over the cell
558: . gradient - The gradient integrated over the cell
559: . cg - A PetscFVCellGeom struct
560: - ctx - A user context
562: Output Parameter:
563: . errInd - The error indicator
565: .seealso: DMAdaptorComputeErrorIndicator()
566: */
567: static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
568: {
569: PetscReal err = 0.;
570: PetscInt c, d;
573: for (c = 0; c < Nc; c++) {
574: for (d = 0; d < dim; ++d) {
575: err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
576: }
577: }
578: *errInd = cg->volume * err;
579: return(0);
580: }
582: static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
583: {
584: PetscDS prob;
585: PetscObject obj;
586: PetscClassId id;
587: void *ctx;
588: PetscQuadrature quad;
589: PetscInt dim, d, cdim, Nc;
590: PetscErrorCode ierr;
593: *errInd = 0.;
594: DMGetDimension(plex, &dim);
595: DMGetCoordinateDim(plex, &cdim);
596: DMGetApplicationContext(plex, &ctx);
597: DMGetDS(plex, &prob);
598: PetscDSGetDiscretization(prob, 0, &obj);
599: PetscObjectGetClassId(obj, &id);
600: if (id == PETSCFV_CLASSID) {
601: const PetscScalar *pointSols;
602: const PetscScalar *pointSol;
603: const PetscScalar *pointGrad;
604: PetscFVCellGeom *cg;
606: PetscFVGetNumComponents((PetscFV) obj, &Nc);
607: VecGetArrayRead(locX, &pointSols);
608: DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);
609: DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);
610: DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);
611: (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);
612: VecRestoreArrayRead(locX, &pointSols);
613: } else {
614: PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
615: PetscFVCellGeom cg;
616: PetscFEGeom fegeom;
617: const PetscReal *quadWeights;
618: PetscReal *coords;
619: PetscInt Nb, fc, Nq, qNc, Nf, f, fieldOffset;
621: fegeom.dim = dim;
622: fegeom.dimEmbed = cdim;
623: PetscDSGetNumFields(prob, &Nf);
624: PetscFEGetQuadrature((PetscFE) obj, &quad);
625: DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);
626: PetscFEGetDimension((PetscFE) obj, &Nb);
627: PetscFEGetNumComponents((PetscFE) obj, &Nc);
628: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
629: PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);
630: PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);
631: DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
632: DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);
633: PetscArrayzero(gradient, cdim*Nc);
634: for (f = 0, fieldOffset = 0; f < Nf; ++f) {
635: PetscInt qc = 0, q;
637: PetscDSGetDiscretization(prob, f, &obj);
638: PetscArrayzero(interpolant,Nc);
639: PetscArrayzero(interpolantGrad, cdim*Nc);
640: for (q = 0; q < Nq; ++q) {
641: PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, x, &fegeom, q, interpolant, interpolantGrad);
642: for (fc = 0; fc < Nc; ++fc) {
643: const PetscReal wt = quadWeights[q*qNc+qc+fc];
645: field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
646: for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
647: }
648: }
649: fieldOffset += Nb;
650: qc += Nc;
651: }
652: PetscFree2(interpolant, interpolantGrad);
653: DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);
654: for (fc = 0; fc < Nc; ++fc) {
655: field[fc] /= cg.volume;
656: for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
657: }
658: (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);
659: PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
660: }
661: return(0);
662: }
664: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
665: {
666: PetscDS prob;
667: void *ctx;
668: MPI_Comm comm;
669: PetscInt numAdapt = adaptor->numSeq, adaptIter;
670: PetscInt dim, coordDim, numFields, cStart, cEnd, cEndInterior, c;
674: DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");
675: VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");
676: PetscObjectGetComm((PetscObject) adaptor, &comm);
677: DMGetDimension(adaptor->idm, &dim);
678: DMGetCoordinateDim(adaptor->idm, &coordDim);
679: DMGetApplicationContext(adaptor->idm, &ctx);
680: DMGetDS(adaptor->idm, &prob);
681: PetscDSGetNumFields(prob, &numFields);
683: /* Adapt until nothing changes */
684: /* Adapt for a specified number of iterates */
685: for (adaptIter = 0; adaptIter < numAdapt-1; ++adaptIter) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));}
686: for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
687: PetscBool adapted = PETSC_FALSE;
688: DM dm = adaptIter ? *adm : adaptor->idm, odm;
689: Vec x = adaptIter ? *ax : inx, locX, ox;
691: DMGetLocalVector(dm, &locX);
692: DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
693: DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);
694: DMAdaptorPreAdapt(adaptor, locX);
695: if (doSolve) {
696: SNES snes;
698: DMAdaptorGetSolver(adaptor, &snes);
699: SNESSolve(snes, NULL, adaptIter ? *ax : x);
700: }
701: /* DMAdaptorMonitor(adaptor);
702: Print iterate, memory used, DM, solution */
703: switch (adaptor->adaptCriterion) {
704: case DM_ADAPTATION_REFINE:
705: DMRefine(dm, comm, &odm);
706: if (!odm) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
707: adapted = PETSC_TRUE;
708: break;
709: case DM_ADAPTATION_LABEL:
710: {
711: /* Adapt DM
712: Create local solution
713: Reconstruct gradients (FVM) or solve adjoint equation (FEM)
714: Produce cellwise error indicator */
715: DM plex;
716: DMLabel adaptLabel;
717: IS refineIS, coarsenIS;
718: Vec errVec;
719: PetscScalar *errArray;
720: const PetscScalar *pointSols;
721: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
722: PetscInt nRefine, nCoarsen;
724: DMConvert(dm, DMPLEX, &plex);
725: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
726: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
727: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
728: cEnd = (cEndInterior < 0) ? cEnd : cEndInterior;
730: VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);
731: VecSetUp(errVec);
732: VecGetArray(errVec, &errArray);
733: VecGetArrayRead(locX, &pointSols);
734: for (c = cStart; c < cEnd; ++c) {
735: PetscReal errInd;
737: DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);
738: errArray[c-cStart] = errInd;
739: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
740: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
741: }
742: VecRestoreArrayRead(locX, &pointSols);
743: VecRestoreArray(errVec, &errArray);
744: PetscGlobalMinMaxReal(PetscObjectComm((PetscObject) adaptor), minMaxInd, minMaxIndGlobal);
745: PetscInfo2(adaptor, "DMAdaptor: error indicator range (%E, %E)\n", minMaxIndGlobal[0], minMaxIndGlobal[1]);
746: /* Compute IS from VecTagger */
747: VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS);
748: VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS);
749: ISGetSize(refineIS, &nRefine);
750: ISGetSize(coarsenIS, &nCoarsen);
751: PetscInfo2(adaptor, "DMAdaptor: numRefine %D, numCoarsen %D\n", nRefine, nCoarsen);
752: if (nRefine) {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS);}
753: if (nCoarsen) {DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);}
754: ISDestroy(&coarsenIS);
755: ISDestroy(&refineIS);
756: VecDestroy(&errVec);
757: /* Adapt DM from label */
758: if (nRefine || nCoarsen) {
759: DMAdaptLabel(dm, adaptLabel, &odm);
760: adapted = PETSC_TRUE;
761: }
762: DMLabelDestroy(&adaptLabel);
763: DMDestroy(&plex);
764: }
765: break;
766: case DM_ADAPTATION_METRIC:
767: {
768: DM dmGrad, dmHess, dmMetric;
769: PetscDS probGrad, probHess;
770: Vec xGrad, xHess, metric;
771: PetscSection sec, msec;
772: PetscScalar *H, *M, integral;
773: PetscReal N;
774: DMLabel bdLabel;
775: PetscInt Nd = coordDim*coordDim, f, vStart, vEnd, v;
777: /* Compute vertexwise gradients from cellwise gradients */
778: DMClone(dm, &dmGrad);
779: DMClone(dm, &dmHess);
780: DMGetDS(dmGrad, &probGrad);
781: DMGetDS(dmHess, &probHess);
782: for (f = 0; f < numFields; ++f) {
783: PetscFE fe, feGrad, feHess;
784: PetscDualSpace Q;
785: DM K;
786: PetscQuadrature q;
787: PetscInt Nc, qorder;
788: const char *prefix;
790: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
791: PetscFEGetNumComponents(fe, &Nc);
792: PetscFEGetDualSpace(fe, &Q);
793: PetscDualSpaceGetDM(Q, &K);
794: DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);
795: PetscFEGetQuadrature(fe, &q);
796: PetscQuadratureGetOrder(q, &qorder);
797: PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);
798: PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feGrad);
799: PetscDSSetDiscretization(probGrad, f, (PetscObject) feGrad);
800: PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feHess);
801: PetscDSSetDiscretization(probHess, f, (PetscObject) feHess);
802: PetscFEDestroy(&feGrad);
803: PetscFEDestroy(&feHess);
804: }
805: DMGetGlobalVector(dmGrad, &xGrad);
806: VecViewFromOptions(x, NULL, "-sol_adapt_loc_pre_view");
807: DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);
808: VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");
809: /* Compute vertexwise Hessians from cellwise Hessians */
810: DMGetGlobalVector(dmHess, &xHess);
811: DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);
812: VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");
813: /* Compute metric */
814: DMClone(dm, &dmMetric);
815: DMGetLocalSection(dm, &sec);
816: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
817: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &msec);
818: PetscSectionSetNumFields(msec, 1);
819: PetscSectionSetFieldComponents(msec, 0, Nd);
820: PetscSectionSetChart(msec, vStart, vEnd);
821: for (v = vStart; v < vEnd; ++v) {
822: PetscSectionSetDof(msec, v, Nd);
823: PetscSectionSetFieldDof(msec, v, 0, Nd);
824: }
825: PetscSectionSetUp(msec);
826: DMSetLocalSection(dmMetric, msec);
827: PetscSectionDestroy(&msec);
828: DMGetLocalVector(dmMetric, &metric);
829: /* N is the target size */
830: N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
831: if (adaptor->monitor) {PetscPrintf(PETSC_COMM_SELF, "N_orig: %D N_adapt: %g\n", vEnd - vStart, N);}
832: /* |H| means take the absolute value of eigenvalues */
833: VecGetArray(xHess, &H);
834: VecGetArray(metric, &M);
835: for (v = vStart; v < vEnd; ++v) {
836: PetscScalar *Hp;
838: DMPlexPointLocalRef(dmHess, v, H, &Hp);
839: DMAdaptorModifyHessian_Private(coordDim, adaptor->h_min, adaptor->h_max, Hp);
840: }
841: /* Pointwise on vertices M(x) = N^{2/d} (\int_\Omega det(|H|)^{p/(2p+d)})^{-2/d} det(|H|)^{-1/(2p+d)} |H| for L_p */
842: PetscDSSetObjective(probHess, 0, detHFunc);
843: DMPlexComputeIntegralFEM(dmHess, xHess, &integral, NULL);
844: for (v = vStart; v < vEnd; ++v) {
845: const PetscInt p = 1;
846: const PetscScalar *Hp;
847: PetscScalar *Mp;
848: PetscReal detH, fact;
849: PetscInt i;
851: DMPlexPointLocalRead(dmHess, v, H, (void *) &Hp);
852: DMPlexPointLocalRef(dmMetric, v, M, &Mp);
853: if (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, Hp);
854: else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, Hp);
855: else SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "Dimension %d not supported", dim);
856: fact = PetscPowReal(N, 2.0/dim) * PetscPowReal(PetscRealPart(integral), -2.0/dim) * PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim));
857: #if 0
858: PetscPrintf(PETSC_COMM_SELF, "fact: %g integral: %g |detH|: %g termA: %g termB: %g\n", fact, integral, PetscAbsReal(detH), PetscPowReal(integral, -2.0/dim), PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim)));
859: DMPrintCellMatrix(v, "H", coordDim, coordDim, Hp);
860: #endif
861: for (i = 0; i < Nd; ++i) {
862: Mp[i] = fact * Hp[i];
863: }
864: }
865: VecRestoreArray(xHess, &H);
866: VecRestoreArray(metric, &M);
867: /* Adapt DM from metric */
868: DMGetLabel(dm, "marker", &bdLabel);
869: DMAdaptMetric(dm, metric, bdLabel, &odm);
870: adapted = PETSC_TRUE;
871: /* Cleanup */
872: DMRestoreLocalVector(dmMetric, &metric);
873: DMDestroy(&dmMetric);
874: DMRestoreGlobalVector(dmHess, &xHess);
875: DMRestoreGlobalVector(dmGrad, &xGrad);
876: DMDestroy(&dmGrad);
877: DMDestroy(&dmHess);
878: }
879: break;
880: default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
881: }
882: DMAdaptorPostAdapt(adaptor);
883: DMRestoreLocalVector(dm, &locX);
884: /* If DM was adapted, replace objects and recreate solution */
885: if (adapted) {
886: const char *name;
888: PetscObjectGetName((PetscObject) dm, &name);
889: PetscObjectSetName((PetscObject) odm, name);
890: /* Reconfigure solver */
891: SNESReset(adaptor->snes);
892: SNESSetDM(adaptor->snes, odm);
893: DMAdaptorSetSolver(adaptor, adaptor->snes);
894: DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);
895: SNESSetFromOptions(adaptor->snes);
896: /* Transfer system */
897: DMCopyDisc(adaptor->idm, odm);
898: /* Transfer solution */
899: DMCreateGlobalVector(odm, &ox);
900: PetscObjectGetName((PetscObject) x, &name);
901: PetscObjectSetName((PetscObject) ox, name);
902: DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);
903: /* Cleanup adaptivity info */
904: if (adaptIter > 0) {PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));}
905: DMForestSetAdaptivityForest(dm, NULL); /* clear internal references to the previous dm */
906: DMDestroy(&dm);
907: VecDestroy(&x);
908: *adm = odm;
909: *ax = ox;
910: } else {
911: *adm = dm;
912: *ax = x;
913: adaptIter = numAdapt;
914: }
915: if (adaptIter < numAdapt-1) {
916: DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");
917: VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");
918: }
919: }
920: DMViewFromOptions(*adm, NULL, "-dm_adapt_view");
921: VecViewFromOptions(*ax, NULL, "-sol_adapt_view");
922: return(0);
923: }
925: /*@
926: DMAdaptorAdapt - Creates a new DM that is adapted to the problem
928: Not collective
930: Input Parameter:
931: + adaptor - The DMAdaptor object
932: . x - The global approximate solution
933: - strategy - The adaptation strategy
935: Output Parameters:
936: + adm - The adapted DM
937: - ax - The adapted solution
939: Options database keys:
940: + -snes_adapt <strategy> : initial, sequential, multigrid
941: . -adapt_gradient_view : View the Clement interpolant of the solution gradient
942: . -adapt_hessian_view : View the Clement interpolant of the solution Hessian
943: - -adapt_metric_view : View the metric tensor for adaptive mesh refinement
945: Note: The available adaptation strategies are:
946: $ 1) Adapt the intial mesh until a quality metric, e,g, a priori error bound, is satisfied
947: $ 2) Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
948: $ 3) Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid
950: Level: intermediate
952: .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
953: @*/
954: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
955: {
959: switch (strategy)
960: {
961: case DM_ADAPTATION_INITIAL:
962: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);
963: break;
964: case DM_ADAPTATION_SEQUENTIAL:
965: DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);
966: break;
967: default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
968: }
969: return(0);
970: }