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