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>
6: #include <petscsnes.h>
7: #include <petscdraw.h>
9: #include <petsc/private/dmadaptorimpl.h>
10: #include <petsc/private/dmpleximpl.h>
11: #include <petsc/private/petscfeimpl.h>
13: PetscClassId DMADAPTOR_CLASSID;
15: PetscFunctionList DMAdaptorList = NULL;
16: PetscBool DMAdaptorRegisterAllCalled = PETSC_FALSE;
18: PetscFunctionList DMAdaptorMonitorList = NULL;
19: PetscFunctionList DMAdaptorMonitorCreateList = NULL;
20: PetscFunctionList DMAdaptorMonitorDestroyList = NULL;
21: PetscBool DMAdaptorMonitorRegisterAllCalled = PETSC_FALSE;
23: const char *const DMAdaptationCriteria[] = {"NONE", "REFINE", "LABEL", "METRIC", "DMAdaptationCriterion", "DM_ADAPTATION_", NULL};
25: /*@C
26: DMAdaptorRegister - Adds a new adaptor component implementation
28: Not Collective
30: Input Parameters:
31: + name - The name of a new user-defined creation routine
32: - create_func - The creation routine
34: Example Usage:
35: .vb
36: DMAdaptorRegister("my_adaptor", MyAdaptorCreate);
37: .ve
39: Then, your adaptor type can be chosen with the procedural interface via
40: .vb
41: DMAdaptorCreate(MPI_Comm, DMAdaptor *);
42: DMAdaptorSetType(DMAdaptor, "my_adaptor");
43: .ve
44: or at runtime via the option
45: .vb
46: -adaptor_type my_adaptor
47: .ve
49: Level: advanced
51: Note:
52: `DMAdaptorRegister()` may be called multiple times to add several user-defined adaptors
54: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorRegisterAll()`, `DMAdaptorRegisterDestroy()`
55: @*/
56: PetscErrorCode DMAdaptorRegister(const char name[], PetscErrorCode (*create_func)(DMAdaptor))
57: {
58: PetscFunctionBegin;
59: PetscCall(DMInitializePackage());
60: PetscCall(PetscFunctionListAdd(&DMAdaptorList, name, create_func));
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor);
65: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor);
67: /*@C
68: DMAdaptorRegisterAll - Registers all of the adaptor components in the `DM` package.
70: Not Collective
72: Level: advanced
74: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorType`, `DMRegisterAll()`, `DMAdaptorRegisterDestroy()`
75: @*/
76: PetscErrorCode DMAdaptorRegisterAll(void)
77: {
78: PetscFunctionBegin;
79: if (DMAdaptorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
80: DMAdaptorRegisterAllCalled = PETSC_TRUE;
82: PetscCall(DMAdaptorRegister(DMADAPTORGRADIENT, DMAdaptorCreate_Gradient));
83: PetscCall(DMAdaptorRegister(DMADAPTORFLUX, DMAdaptorCreate_Flux));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /*@C
88: DMAdaptorRegisterDestroy - This function destroys the registered `DMAdaptorType`. It is called from `PetscFinalize()`.
90: Not collective
92: Level: developer
94: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorRegisterAll()`, `DMAdaptorType`, `PetscFinalize()`
95: @*/
96: PetscErrorCode DMAdaptorRegisterDestroy(void)
97: {
98: PetscFunctionBegin;
99: PetscCall(PetscFunctionListDestroy(&DMAdaptorList));
100: DMAdaptorRegisterAllCalled = PETSC_FALSE;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode DMAdaptorMonitorMakeKey_Internal(const char name[], PetscViewerType vtype, PetscViewerFormat format, char key[])
105: {
106: PetscFunctionBegin;
107: PetscCall(PetscStrncpy(key, name, PETSC_MAX_PATH_LEN));
108: PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
109: PetscCall(PetscStrlcat(key, vtype, PETSC_MAX_PATH_LEN));
110: PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
111: PetscCall(PetscStrlcat(key, PetscViewerFormats[format], PETSC_MAX_PATH_LEN));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: /*@C
116: DMAdaptorMonitorRegister - Registers a mesh adaptation monitor routine that may be accessed with `DMAdaptorMonitorSetFromOptions()`
118: Not Collective
120: Input Parameters:
121: + name - name of a new monitor routine
122: . vtype - A `PetscViewerType` for the output
123: . format - A `PetscViewerFormat` for the output
124: . monitor - Monitor routine
125: . create - Creation routine, or `NULL`
126: - destroy - Destruction routine, or `NULL`
128: Level: advanced
130: Note:
131: `DMAdaptorMonitorRegister()` may be called multiple times to add several user-defined monitors.
133: Example Usage:
134: .vb
135: DMAdaptorMonitorRegister("my_monitor", PETSCVIEWERASCII, PETSC_VIEWER_ASCII_INFO_DETAIL, MyMonitor, NULL, NULL);
136: .ve
138: Then, your monitor can be chosen with the procedural interface via
139: .vb
140: DMAdaptorMonitorSetFromOptions(ksp, "-adaptor_monitor_my_monitor", "my_monitor", NULL)
141: .ve
142: or at runtime via the option `-adaptor_monitor_my_monitor`
144: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorRegisterAll()`, `DMAdaptorMonitorSetFromOptions()`
145: @*/
146: PetscErrorCode DMAdaptorMonitorRegister(const char name[], PetscViewerType vtype, PetscViewerFormat format, PetscErrorCode (*monitor)(DMAdaptor, PetscInt, DM, DM, PetscInt, PetscReal[], Vec, PetscViewerAndFormat *), PetscErrorCode (*create)(PetscViewer, PetscViewerFormat, void *, PetscViewerAndFormat **), PetscErrorCode (*destroy)(PetscViewerAndFormat **))
147: {
148: char key[PETSC_MAX_PATH_LEN];
150: PetscFunctionBegin;
151: PetscCall(SNESInitializePackage());
152: PetscCall(DMAdaptorMonitorMakeKey_Internal(name, vtype, format, key));
153: PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorList, key, monitor));
154: if (create) PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorCreateList, key, create));
155: if (destroy) PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorDestroyList, key, destroy));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@C
160: DMAdaptorMonitorRegisterDestroy - This function destroys the registered monitors for `DMAdaptor`. It is called from `PetscFinalize()`.
162: Not collective
164: Level: developer
166: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorMonitorRegisterAll()`, `DMAdaptor`, `PetscFinalize()`
167: @*/
168: PetscErrorCode DMAdaptorMonitorRegisterDestroy(void)
169: {
170: PetscFunctionBegin;
171: PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorList));
172: PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorCreateList));
173: PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorDestroyList));
174: DMAdaptorMonitorRegisterAllCalled = PETSC_FALSE;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: DMAdaptorCreate - Create a `DMAdaptor` object. Its purpose is to construct a adaptation `DMLabel` or metric `Vec` that can be used to modify the `DM`.
181: Collective
183: Input Parameter:
184: . comm - The communicator for the `DMAdaptor` object
186: Output Parameter:
187: . adaptor - The `DMAdaptor` object
189: Level: beginner
191: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`, `PetscConvEst`, `PetscConvEstCreate()`
192: @*/
193: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
194: {
195: VecTaggerBox refineBox, coarsenBox;
197: PetscFunctionBegin;
198: PetscAssertPointer(adaptor, 2);
199: PetscCall(PetscSysInitializePackage());
201: PetscCall(PetscHeaderCreate(*adaptor, DMADAPTOR_CLASSID, "DMAdaptor", "DM Adaptor", "DMAdaptor", comm, DMAdaptorDestroy, DMAdaptorView));
202: (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
203: (*adaptor)->numSeq = 1;
204: (*adaptor)->Nadapt = -1;
205: (*adaptor)->refinementFactor = 2.0;
206: refineBox.min = refineBox.max = PETSC_MAX_REAL;
207: PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->refineTag));
208: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->refineTag, "refine_"));
209: PetscCall(VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE));
210: PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox));
211: coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
212: PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->coarsenTag));
213: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->coarsenTag, "coarsen_"));
214: PetscCall(VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE));
215: PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: DMAdaptorDestroy - Destroys a `DMAdaptor` object
222: Collective
224: Input Parameter:
225: . adaptor - The `DMAdaptor` object
227: Level: beginner
229: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
230: @*/
231: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
232: {
233: PetscFunctionBegin;
234: if (!*adaptor) PetscFunctionReturn(PETSC_SUCCESS);
236: if (--((PetscObject)*adaptor)->refct > 0) {
237: *adaptor = NULL;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
240: PetscCall(VecTaggerDestroy(&(*adaptor)->refineTag));
241: PetscCall(VecTaggerDestroy(&(*adaptor)->coarsenTag));
242: PetscCall(PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx));
243: PetscCall(DMAdaptorMonitorCancel(*adaptor));
244: PetscCall(PetscHeaderDestroy(adaptor));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: DMAdaptorSetType - Sets the particular implementation for a adaptor.
251: Collective
253: Input Parameters:
254: + adaptor - The `DMAdaptor`
255: - method - The name of the adaptor type
257: Options Database Key:
258: . -adaptor_type type - Sets the adaptor type; see `DMAdaptorType`
260: Level: intermediate
262: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorGetType()`, `DMAdaptorCreate()`
263: @*/
264: PetscErrorCode DMAdaptorSetType(DMAdaptor adaptor, DMAdaptorType method)
265: {
266: PetscErrorCode (*r)(DMAdaptor);
267: PetscBool match;
269: PetscFunctionBegin;
271: PetscCall(PetscObjectTypeCompare((PetscObject)adaptor, method, &match));
272: if (match) PetscFunctionReturn(PETSC_SUCCESS);
274: PetscCall(DMAdaptorRegisterAll());
275: PetscCall(PetscFunctionListFind(DMAdaptorList, method, &r));
276: PetscCheck(r, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMAdaptor type: %s", method);
278: PetscTryTypeMethod(adaptor, destroy);
279: PetscCall(PetscMemzero(adaptor->ops, sizeof(*adaptor->ops)));
280: PetscCall(PetscObjectChangeTypeName((PetscObject)adaptor, method));
281: PetscCall((*r)(adaptor));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: DMAdaptorGetType - Gets the type name (as a string) from the adaptor.
288: Not Collective
290: Input Parameter:
291: . adaptor - The `DMAdaptor`
293: Output Parameter:
294: . type - The `DMAdaptorType` name
296: Level: intermediate
298: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorSetType()`, `DMAdaptorCreate()`
299: @*/
300: PetscErrorCode DMAdaptorGetType(DMAdaptor adaptor, DMAdaptorType *type)
301: {
302: PetscFunctionBegin;
304: PetscAssertPointer(type, 2);
305: PetscCall(DMAdaptorRegisterAll());
306: *type = ((PetscObject)adaptor)->type_name;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode PetscViewerAndFormatCreate_Internal(PetscViewer viewer, PetscViewerFormat format, PetscCtx ctx, PetscViewerAndFormat **vf)
311: {
312: PetscFunctionBegin;
313: PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
314: (*vf)->data = ctx;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@C
319: DMAdaptorMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
320: the error etc.
322: Logically Collective
324: Input Parameters:
325: + adaptor - the `DMAdaptor`
326: . monitor - pointer to function (if this is `NULL`, it turns off monitoring
327: . ctx - [optional] context for private data for the monitor routine (use `NULL` if no context is needed)
328: - monitordestroy - [optional] routine that frees monitor context (may be `NULL`), see `PetscCtxDestroyFn` for its calling sequence
330: Calling sequence of `monitor`:
331: + adaptor - the `DMAdaptor`
332: . it - iteration number
333: . odm - the original `DM`
334: . adm - the adapted `DM`
335: . Nf - number of fields
336: . enorms - (estimated) 2-norm of the error for each field
337: . error - `Vec` of cellwise errors
338: - ctx - optional monitoring context, as set by `DMAdaptorMonitorSet()`
340: Options Database Keys:
341: + -adaptor_monitor_size - sets `DMAdaptorMonitorSize()`
342: . -adaptor_monitor_error - sets `DMAdaptorMonitorError()`
343: . -adaptor_monitor_error draw - sets `DMAdaptorMonitorErrorDraw()` and plots error
344: . -adaptor_monitor_error draw::draw_lg - sets `DMAdaptorMonitorErrorDrawLG()` and plots error
345: - -dm_adaptor_monitor_cancel - Cancels all monitors that have been hardwired into a code by calls to `DMAdaptorMonitorSet()`, but does not cancel those set via the options database.
347: Level: beginner
349: .seealso: [](ch_snes), `DMAdaptorMonitorError()`, `DMAdaptor`, `PetscCtxDestroyFn`
350: @*/
351: PetscErrorCode DMAdaptorMonitorSet(DMAdaptor adaptor, PetscErrorCode (*monitor)(DMAdaptor adaptor, PetscInt it, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscCtx ctx), PetscCtx ctx, PetscCtxDestroyFn *monitordestroy)
352: {
353: PetscFunctionBegin;
355: for (PetscInt i = 0; i < adaptor->numbermonitors; i++) {
356: PetscBool identical;
358: PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))(PetscVoidFn *)monitor, ctx, monitordestroy, (PetscErrorCode (*)(void))(PetscVoidFn *)adaptor->monitor[i], adaptor->monitorcontext[i], adaptor->monitordestroy[i], &identical));
359: if (identical) PetscFunctionReturn(PETSC_SUCCESS);
360: }
361: PetscCheck(adaptor->numbermonitors < MAXDMADAPTORMONITORS, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_OUTOFRANGE, "Too many DMAdaptor monitors set");
362: adaptor->monitor[adaptor->numbermonitors] = monitor;
363: adaptor->monitordestroy[adaptor->numbermonitors] = monitordestroy;
364: adaptor->monitorcontext[adaptor->numbermonitors++] = ctx;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@
369: DMAdaptorMonitorCancel - Clears all monitors for a `DMAdaptor` object.
371: Logically Collective
373: Input Parameter:
374: . adaptor - the `DMAdaptor`
376: Options Database Key:
377: . -dm_adaptor_monitor_cancel - Cancels all monitors that have been hardwired into a code by calls to `DMAdaptorMonitorSet()`, but does not cancel those set via the options database.
379: Level: intermediate
381: .seealso: [](ch_snes), `DMAdaptorMonitorError()`, `DMAdaptorMonitorSet()`, `DMAdaptor`
382: @*/
383: PetscErrorCode DMAdaptorMonitorCancel(DMAdaptor adaptor)
384: {
385: PetscFunctionBegin;
387: for (PetscInt i = 0; i < adaptor->numbermonitors; ++i) {
388: if (adaptor->monitordestroy[i]) PetscCall((*adaptor->monitordestroy[i])(&adaptor->monitorcontext[i]));
389: }
390: adaptor->numbermonitors = 0;
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@C
395: DMAdaptorMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user in the options database
397: Collective
399: Input Parameters:
400: + adaptor - `DMadaptor` object you wish to monitor
401: . opt - the command line option for this monitor
402: . name - the monitor type one is seeking
403: - ctx - An optional application context for the monitor, or `NULL`
405: Level: developer
407: .seealso: [](ch_snes), `DMAdaptorMonitorRegister()`, `DMAdaptorMonitorSet()`, `PetscOptionsGetViewer()`
408: @*/
409: PetscErrorCode DMAdaptorMonitorSetFromOptions(DMAdaptor adaptor, const char opt[], const char name[], PetscCtx ctx)
410: {
411: PetscErrorCode (*mfunc)(DMAdaptor, PetscInt, DM, DM, PetscInt, PetscReal[], Vec, void *);
412: PetscErrorCode (*cfunc)(PetscViewer, PetscViewerFormat, void *, PetscViewerAndFormat **);
413: PetscErrorCode (*dfunc)(PetscViewerAndFormat **);
414: PetscViewerAndFormat *vf;
415: PetscViewer viewer;
416: PetscViewerFormat format;
417: PetscViewerType vtype;
418: char key[PETSC_MAX_PATH_LEN];
419: PetscBool flg;
420: const char *prefix = NULL;
422: PetscFunctionBegin;
423: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
424: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)adaptor), ((PetscObject)adaptor)->options, prefix, opt, &viewer, &format, &flg));
425: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
427: PetscCall(PetscViewerGetType(viewer, &vtype));
428: PetscCall(DMAdaptorMonitorMakeKey_Internal(name, vtype, format, key));
429: PetscCall(PetscFunctionListFind(DMAdaptorMonitorList, key, &mfunc));
430: PetscCall(PetscFunctionListFind(DMAdaptorMonitorCreateList, key, &cfunc));
431: PetscCall(PetscFunctionListFind(DMAdaptorMonitorDestroyList, key, &dfunc));
432: if (!cfunc) cfunc = PetscViewerAndFormatCreate_Internal;
433: if (!dfunc) dfunc = PetscViewerAndFormatDestroy;
435: PetscCall((*cfunc)(viewer, format, ctx, &vf));
436: PetscCall(PetscViewerDestroy(&viewer));
437: PetscCall(DMAdaptorMonitorSet(adaptor, mfunc, vf, (PetscCtxDestroyFn *)dfunc));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: /*@
442: DMAdaptorSetOptionsPrefix - Sets the prefix used for searching for all `DMAdaptor` options in the database.
444: Logically Collective
446: Input Parameters:
447: + adaptor - the `DMAdaptor`
448: - prefix - the prefix to prepend to all option names
450: Level: advanced
452: Note:
453: A hyphen (-) must NOT be given at the beginning of the prefix name.
454: The first character of all runtime options is AUTOMATICALLY the hyphen.
456: .seealso: [](ch_snes), `DMAdaptor`, `SNESSetOptionsPrefix()`, `DMAdaptorSetFromOptions()`
457: @*/
458: PetscErrorCode DMAdaptorSetOptionsPrefix(DMAdaptor adaptor, const char prefix[])
459: {
460: PetscFunctionBegin;
462: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor, prefix));
463: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor->refineTag, prefix));
464: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)adaptor->refineTag, "refine_"));
465: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor->coarsenTag, prefix));
466: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)adaptor->coarsenTag, "coarsen_"));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: DMAdaptorSetFromOptions - Sets properties of a `DMAdaptor` object from values in the options database
473: Collective
475: Input Parameter:
476: . adaptor - The `DMAdaptor` object
478: Options Database Keys:
479: + -adaptor_monitor_size - Monitor the mesh size
480: . -adaptor_monitor_error - Monitor the solution error
481: . -adaptor_sequence_num num - Number of adaptations to generate an optimal grid
482: . -adaptor_target_num num - Set the target number of vertices N_adapt, -1 for automatic determination
483: . -adaptor_refinement_factor r - Set r such that N_adapt = r^dim N_orig
484: - -adaptor_mixed_setup_function func - Set the function func that sets up the mixed problem
486: Level: beginner
488: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
489: @*/
490: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
491: {
492: char typeName[PETSC_MAX_PATH_LEN];
493: const char *defName = DMADAPTORGRADIENT;
494: char funcname[PETSC_MAX_PATH_LEN];
495: DMAdaptationCriterion criterion = DM_ADAPTATION_NONE;
496: PetscBool flg;
498: PetscFunctionBegin;
499: PetscObjectOptionsBegin((PetscObject)adaptor);
500: PetscCall(PetscOptionsFList("-adaptor_type", "DMAdaptor", "DMAdaptorSetType", DMAdaptorList, defName, typeName, 1024, &flg));
501: if (flg) PetscCall(DMAdaptorSetType(adaptor, typeName));
502: else if (!((PetscObject)adaptor)->type_name) PetscCall(DMAdaptorSetType(adaptor, defName));
503: PetscCall(PetscOptionsEnum("-adaptor_criterion", "Criterion used to drive adaptation", "", DMAdaptationCriteria, (PetscEnum)criterion, (PetscEnum *)&criterion, &flg));
504: if (flg) PetscCall(DMAdaptorSetCriterion(adaptor, criterion));
505: PetscCall(PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL));
506: PetscCall(PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL));
507: PetscCall(PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL));
508: PetscCall(PetscOptionsString("-adaptor_mixed_setup_function", "Function to setup the mixed problem", "DMAdaptorSetMixedSetupFunction", funcname, funcname, sizeof(funcname), &flg));
509: if (flg) {
510: PetscErrorCode (*setupFunc)(DMAdaptor, DM);
512: PetscCall(PetscDLSym(NULL, funcname, (void **)&setupFunc));
513: PetscCheck(setupFunc, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
514: PetscCall(DMAdaptorSetMixedSetupFunction(adaptor, setupFunc));
515: }
516: PetscCall(DMAdaptorMonitorSetFromOptions(adaptor, "-adaptor_monitor_size", "size", adaptor));
517: PetscCall(DMAdaptorMonitorSetFromOptions(adaptor, "-adaptor_monitor_error", "error", adaptor));
518: PetscOptionsEnd();
519: PetscCall(VecTaggerSetFromOptions(adaptor->refineTag));
520: PetscCall(VecTaggerSetFromOptions(adaptor->coarsenTag));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: DMAdaptorView - Views a `DMAdaptor` object
527: Collective
529: Input Parameters:
530: + adaptor - The `DMAdaptor` object
531: - viewer - The `PetscViewer` object
533: Level: beginner
535: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
536: @*/
537: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
538: {
539: PetscFunctionBegin;
540: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer));
541: PetscCall(PetscViewerASCIIPrintf(viewer, "DM Adaptor\n"));
542: PetscCall(PetscViewerASCIIPrintf(viewer, " sequence length: %" PetscInt_FMT "\n", adaptor->numSeq));
543: PetscCall(VecTaggerView(adaptor->refineTag, viewer));
544: PetscCall(VecTaggerView(adaptor->coarsenTag, viewer));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: /*@
549: DMAdaptorGetSolver - Gets the solver used to produce discrete solutions
551: Not Collective
553: Input Parameter:
554: . adaptor - The `DMAdaptor` object
556: Output Parameter:
557: . snes - The solver
559: Level: intermediate
561: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
562: @*/
563: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
564: {
565: PetscFunctionBegin;
567: PetscAssertPointer(snes, 2);
568: *snes = adaptor->snes;
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@
573: DMAdaptorSetSolver - Sets the solver used to produce discrete solutions
575: Not Collective
577: Input Parameters:
578: + adaptor - The `DMAdaptor` object
579: - snes - The solver, this MUST have an attached `DM`/`PetscDS`, so that the exact solution can be computed
581: Level: intermediate
583: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
584: @*/
585: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
586: {
587: PetscFunctionBegin;
590: adaptor->snes = snes;
591: PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: /*@
596: DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter
598: Not Collective
600: Input Parameter:
601: . adaptor - The `DMAdaptor` object
603: Output Parameter:
604: . num - The number of adaptations
606: Level: intermediate
608: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
609: @*/
610: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
611: {
612: PetscFunctionBegin;
614: PetscAssertPointer(num, 2);
615: *num = adaptor->numSeq;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@
620: DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
622: Not Collective
624: Input Parameters:
625: + adaptor - The `DMAdaptor` object
626: - num - The number of adaptations
628: Level: intermediate
630: .seealso: [](ch_dmbase), `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
631: @*/
632: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
633: {
634: PetscFunctionBegin;
636: adaptor->numSeq = num;
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, PetscCtx ctx)
641: {
642: PetscFunctionBeginUser;
643: PetscCall(DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: /*@
648: DMAdaptorSetUp - After the solver is specified, creates data structures for controlling adaptivity
650: Collective
652: Input Parameter:
653: . adaptor - The `DMAdaptor` object
655: Level: beginner
657: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
658: @*/
659: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
660: {
661: PetscDS prob;
662: PetscInt Nf;
664: PetscFunctionBegin;
665: PetscCall(VecTaggerSetUp(adaptor->refineTag));
666: PetscCall(VecTaggerSetUp(adaptor->coarsenTag));
667: PetscCall(DMGetDS(adaptor->idm, &prob));
668: PetscCall(PetscDSGetNumFields(prob, &Nf));
669: PetscCall(PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx));
670: for (PetscInt f = 0; f < Nf; ++f) {
671: PetscCall(PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]));
672: /* TODO Have a flag that forces projection rather than using the exact solution */
673: if (adaptor->exactSol[0]) PetscCall(DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private));
674: }
675: PetscTryTypeMethod(adaptor, setup);
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }
679: /*@C
680: DMAdaptorGetTransferFunction - Get the callback used by a `DMAdaptor` to transfer a solution vector from an old `DM` to the adapted `DM`
682: Not Collective
684: Input Parameter:
685: . adaptor - the `DMAdaptor` object
687: Output Parameter:
688: . tfunc - pointer to the transfer callback
690: Calling sequence of `tfunc`:
691: + adaptor - the `DMAdaptor` object
692: . dm - the current `DM`
693: . xin - the current solution
694: . newdm - the adapted `DM`
695: . xout - the transferred solution on `newdm`
696: - ctx - application context, set with `DMSetApplicationContext()`
698: Level: developer
700: .seealso: `DMAdaptor`, `DMAdaptorSetTransferFunction()`, `DMAdaptorAdapt()`
701: @*/
702: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor adaptor, DM dm, Vec xin, DM newdm, Vec xout, PetscCtx ctx))
703: {
704: PetscFunctionBegin;
705: *tfunc = adaptor->ops->transfersolution;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /*@C
710: DMAdaptorSetTransferFunction - Set the callback used by a `DMAdaptor` to transfer a solution vector from an old `DM` to the adapted `DM`
712: Logically Collective
714: Input Parameters:
715: + adaptor - the `DMAdaptor` object
716: - tfunc - the transfer callback
718: Calling sequence of `tfunc`:
719: + adaptor - the `DMAdaptor` object
720: . dm - the current `DM`
721: . xin - the current solution
722: . newdm - the adapted `DM`
723: . xout - the transferred solution on `newdm`
724: - ctx - application context, set with `DMSetApplicationContext()`
726: Level: developer
728: .seealso: `DMAdaptor`, `DMAdaptorGetTransferFunction()`, `DMAdaptorAdapt()`
729: @*/
730: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor adaptor, DM dm, Vec xin, DM newdm, Vec xout, PetscCtx ctx))
731: {
732: PetscFunctionBegin;
733: adaptor->ops->transfersolution = tfunc;
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: static PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
738: {
739: DM plex;
740: PetscDS prob;
741: PetscObject obj;
742: PetscClassId id;
743: PetscBool isForest;
745: PetscFunctionBegin;
746: PetscCall(DMConvert(adaptor->idm, DMPLEX, &plex));
747: PetscCall(DMGetDS(adaptor->idm, &prob));
748: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
749: PetscCall(PetscObjectGetClassId(obj, &id));
750: PetscCall(DMIsForest(adaptor->idm, &isForest));
751: if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
752: if (isForest) adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
753: #if PetscDefined(HAVE_PRAGMATIC)
754: else {
755: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
756: }
757: #elif PetscDefined(HAVE_MMG)
758: else {
759: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
760: }
761: #elif PetscDefined(HAVE_PARMMG)
762: else {
763: adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
764: }
765: #else
766: else {
767: adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
768: }
769: #endif
770: }
771: if (id == PETSCFV_CLASSID) {
772: adaptor->femType = PETSC_FALSE;
773: } else {
774: adaptor->femType = PETSC_TRUE;
775: }
776: if (adaptor->femType) {
777: /* Compute local solution bc */
778: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
779: } else {
780: PetscFV fvm = (PetscFV)obj;
781: PetscLimiter noneLimiter;
782: Vec grad;
784: PetscCall(PetscFVGetComputeGradients(fvm, &adaptor->computeGradient));
785: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
786: /* Use no limiting when reconstructing gradients for adaptivity */
787: PetscCall(PetscFVGetLimiter(fvm, &adaptor->limiter));
788: PetscCall(PetscObjectReference((PetscObject)adaptor->limiter));
789: PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
790: PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
791: PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
792: /* Get FVM data */
793: PetscCall(DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM));
794: PetscCall(VecGetDM(adaptor->cellGeom, &adaptor->cellDM));
795: PetscCall(VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
796: /* Compute local solution bc */
797: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
798: /* Compute gradients */
799: PetscCall(DMCreateGlobalVector(adaptor->gradDM, &grad));
800: PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
801: PetscCall(DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad));
802: PetscCall(DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
803: PetscCall(DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
804: PetscCall(VecDestroy(&grad));
805: PetscCall(VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
806: }
807: PetscCall(DMDestroy(&plex));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
812: {
813: PetscReal time = 0.0;
814: Mat interp;
815: void *ctx;
817: PetscFunctionBegin;
818: PetscCall(DMGetApplicationContext(dm, &ctx));
819: if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
820: else {
821: switch (adaptor->adaptCriterion) {
822: case DM_ADAPTATION_LABEL:
823: PetscCall(DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time));
824: break;
825: case DM_ADAPTATION_REFINE:
826: case DM_ADAPTATION_METRIC:
827: PetscCall(DMCreateInterpolation(dm, adm, &interp, NULL));
828: PetscCall(MatInterpolate(interp, x, ax));
829: PetscCall(DMInterpolate(dm, interp, adm));
830: PetscCall(MatDestroy(&interp));
831: break;
832: default:
833: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
834: }
835: }
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
840: {
841: PetscDS prob;
842: PetscObject obj;
843: PetscClassId id;
845: PetscFunctionBegin;
846: PetscCall(DMGetDS(adaptor->idm, &prob));
847: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
848: PetscCall(PetscObjectGetClassId(obj, &id));
849: if (id == PETSCFV_CLASSID) {
850: PetscFV fvm = (PetscFV)obj;
852: PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
853: /* Restore original limiter */
854: PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));
856: PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
857: PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
858: PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
859: }
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: /*
864: DMAdaptorComputeCellErrorIndicator_Gradient - Use the integrated gradient as an error indicator in the `DMAdaptor`
866: Input Parameters:
867: + adaptor - The `DMAdaptor` object
868: . dim - The topological dimension
869: . cell - The cell
870: . field - The field integrated over the cell
871: . gradient - The gradient integrated over the cell
872: . cg - A `PetscFVCellGeom` struct
873: - ctx - An application context
875: Output Parameter:
876: . errInd - The error indicator
878: Developer Note:
879: Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API
881: .seealso: [](ch_dmbase), `DMAdaptor`
882: */
883: static PetscErrorCode DMAdaptorComputeCellErrorIndicator_Gradient(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, PetscCtx ctx)
884: {
885: PetscReal err = 0.;
887: PetscFunctionBeginHot;
888: for (PetscInt c = 0; c < Nc; c++) {
889: for (PetscInt d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
890: }
891: *errInd = cg->volume * err;
892: PetscFunctionReturn(PETSC_SUCCESS);
893: }
895: static PetscErrorCode DMAdaptorComputeErrorIndicator_Gradient(DMAdaptor adaptor, Vec locX, Vec errVec)
896: {
897: DM dm, plex, edm, eplex;
898: PetscDS ds;
899: PetscObject obj;
900: PetscClassId id;
901: void *ctx;
902: PetscQuadrature quad;
903: PetscScalar *earray;
904: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL};
905: PetscInt dim, cdim, cStart, cEnd, Nf, Nc;
907: PetscFunctionBegin;
908: PetscCall(VecGetDM(locX, &dm));
909: PetscCall(DMConvert(dm, DMPLEX, &plex));
910: PetscCall(VecGetDM(errVec, &edm));
911: PetscCall(DMConvert(edm, DMPLEX, &eplex));
912: PetscCall(DMGetDimension(plex, &dim));
913: PetscCall(DMGetCoordinateDim(plex, &cdim));
914: PetscCall(DMGetApplicationContext(plex, &ctx));
915: PetscCall(DMGetDS(plex, &ds));
916: PetscCall(PetscDSGetNumFields(ds, &Nf));
917: PetscCall(PetscDSGetDiscretization(ds, 0, &obj));
918: PetscCall(PetscObjectGetClassId(obj, &id));
920: PetscCall(VecGetArray(errVec, &earray));
921: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
922: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
923: PetscScalar *eval;
924: PetscReal errInd = 0.;
926: if (id == PETSCFV_CLASSID) {
927: PetscFV fv = (PetscFV)obj;
928: const PetscScalar *pointSols;
929: const PetscScalar *pointSol;
930: const PetscScalar *pointGrad;
931: PetscFVCellGeom *cg;
933: PetscCall(PetscFVGetNumComponents(fv, &Nc));
934: PetscCall(VecGetArrayRead(locX, &pointSols));
935: PetscCall(DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol));
936: PetscCall(DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad));
937: PetscCall(DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg));
938: PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, pointSol, pointGrad, cg, &errInd, ctx);
939: PetscCall(VecRestoreArrayRead(locX, &pointSols));
940: } else {
941: PetscFE fe = (PetscFE)obj;
942: PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
943: PetscFVCellGeom cg;
944: PetscFEGeom fegeom;
945: const PetscReal *quadWeights;
946: PetscReal *coords;
947: PetscInt Nb, Nq, qNc;
949: fegeom.dim = dim;
950: fegeom.dimEmbed = cdim;
951: PetscCall(PetscFEGetNumComponents(fe, &Nc));
952: PetscCall(PetscFEGetQuadrature(fe, &quad));
953: PetscCall(PetscFEGetDimension(fe, &Nb));
954: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights));
955: PetscCall(PetscCalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ));
956: PetscCall(PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad));
957: PetscCall(DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
958: PetscCall(DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL));
959: PetscCall(PetscArrayzero(gradient, cdim * Nc));
960: PetscCall(DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x));
961: for (PetscInt f = 0; f < Nf; ++f) {
962: PetscInt qc = 0;
964: PetscCall(PetscDSGetDiscretization(ds, f, &obj));
965: PetscCall(PetscArrayzero(interpolant, Nc));
966: PetscCall(PetscArrayzero(interpolantGrad, cdim * Nc));
967: for (PetscInt q = 0; q < Nq; ++q) {
968: PetscCall(PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad));
969: for (PetscInt fc = 0; fc < Nc; ++fc) {
970: const PetscReal wt = quadWeights[q * qNc + qc + fc];
972: field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
973: for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
974: }
975: }
976: qc += Nc;
977: }
978: PetscCall(PetscFree2(interpolant, interpolantGrad));
979: PetscCall(DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x));
980: for (PetscInt fc = 0; fc < Nc; ++fc) {
981: field[fc] /= cg.volume;
982: for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
983: }
984: PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, field, gradient, &cg, &errInd, ctx);
985: PetscCall(PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
986: }
987: PetscCall(DMPlexPointGlobalRef(eplex, cell, earray, (void *)&eval));
988: eval[0] = errInd;
989: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
990: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
991: }
992: PetscCall(VecRestoreArray(errVec, &earray));
993: PetscCall(DMDestroy(&plex));
994: PetscCall(DMDestroy(&eplex));
995: PetscCall(PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxInd));
996: PetscCall(PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxInd[0], (double)minMaxInd[1]));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: static PetscErrorCode DMAdaptorComputeErrorIndicator_Flux(DMAdaptor adaptor, Vec lu, Vec errVec)
1001: {
1002: DM dm, mdm;
1003: SNES msnes;
1004: Vec mu, lmu;
1005: void *ctx;
1006: const char *prefix;
1008: PetscFunctionBegin;
1009: PetscCall(VecGetDM(lu, &dm));
1011: // Set up and solve mixed problem
1012: PetscCall(DMClone(dm, &mdm));
1013: PetscCall(SNESCreate(PetscObjectComm((PetscObject)mdm), &msnes));
1014: PetscCall(SNESSetDM(msnes, mdm));
1015: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
1016: PetscCall(SNESSetOptionsPrefix(msnes, prefix));
1017: PetscCall(SNESAppendOptionsPrefix(msnes, "mixed_"));
1019: PetscTryTypeMethod(adaptor, mixedsetup, mdm);
1020: PetscCall(DMGetApplicationContext(dm, &ctx));
1021: PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, ctx));
1022: PetscCall(SNESSetFromOptions(msnes));
1024: PetscCall(DMCreateGlobalVector(mdm, &mu));
1025: PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution"));
1026: PetscCall(SNESSolve(msnes, NULL, mu));
1027: PetscCall(VecViewFromOptions(mu, (PetscObject)adaptor, "-adapt_mixed_sol_vec_view"));
1029: PetscCall(DMGetLocalVector(mdm, &lmu));
1030: PetscCall(DMGlobalToLocal(mdm, mu, INSERT_VALUES, lmu));
1031: PetscCall(DMPlexInsertBoundaryValues(mdm, PETSC_TRUE, lmu, 0.0, NULL, NULL, NULL));
1032: PetscCall(DMPlexComputeL2FluxDiffVecLocal(lu, 0, lmu, 0, errVec));
1033: PetscCall(DMRestoreLocalVector(mdm, &lmu));
1034: PetscCall(VecDestroy(&mu));
1035: PetscCall(SNESDestroy(&msnes));
1036: PetscCall(DMDestroy(&mdm));
1037: PetscFunctionReturn(PETSC_SUCCESS);
1038: }
1040: /*@
1041: DMAdaptorMonitor - runs the user provided monitor routines, if they exist
1043: Collective
1045: Input Parameters:
1046: + adaptor - the `DMAdaptor`
1047: . it - iteration number
1048: . odm - the original `DM`
1049: . adm - the adapted `DM`
1050: . Nf - the number of fields
1051: . enorms - the 2-norm error values for each field
1052: - error - `Vec` of cellwise errors
1054: Level: developer
1056: Note:
1057: This routine is called by the `DMAdaptor` implementations.
1058: It does not typically need to be called by the user.
1060: .seealso: [](ch_snes), `DMAdaptorMonitorSet()`
1061: @*/
1062: PetscErrorCode DMAdaptorMonitor(DMAdaptor adaptor, PetscInt it, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error)
1063: {
1064: PetscFunctionBegin;
1065: for (PetscInt i = 0; i < adaptor->numbermonitors; ++i) PetscCall((*adaptor->monitor[i])(adaptor, it, odm, adm, Nf, enorms, error, adaptor->monitorcontext[i]));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@C
1070: DMAdaptorMonitorSize - Prints the mesh sizes at each iteration of an adaptation loop.
1072: Collective
1074: Input Parameters:
1075: + adaptor - the `DMAdaptor`
1076: . n - iteration number
1077: . odm - the original `DM`
1078: . adm - the adapted `DM`
1079: . Nf - number of fields
1080: . enorms - 2-norm error values for each field (may be estimated).
1081: . error - `Vec` of cellwise errors
1082: - vf - The viewer context
1084: Options Database Key:
1085: . -adaptor_monitor_size - Activates `DMAdaptorMonitorSize()`
1087: Level: intermediate
1089: Note:
1090: This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
1091: to be used during the adaptation loop.
1093: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorError()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorErrorDrawLG()`
1094: @*/
1095: PetscErrorCode DMAdaptorMonitorSize(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
1096: {
1097: PetscViewer viewer = vf->viewer;
1098: PetscViewerFormat format = vf->format;
1099: PetscInt tablevel, cStart, cEnd, acStart, acEnd;
1100: const char *prefix;
1101: PetscMPIInt rank;
1103: PetscFunctionBegin;
1105: PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
1106: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
1107: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)adaptor), &rank));
1108: PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
1109: PetscCall(DMPlexGetHeightStratum(adm, 0, &acStart, &acEnd));
1111: PetscCall(PetscViewerPushFormat(viewer, format));
1112: PetscCall(PetscViewerASCIIAddTab(viewer, tablevel));
1113: if (n == 0 && prefix) PetscCall(PetscViewerASCIIPrintf(viewer, " Sizes for %s adaptation.\n", prefix));
1114: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " DMAdaptor rank %d N_orig: %" PetscInt_FMT " N_adapt: %" PetscInt_FMT "\n", n, rank, cEnd - cStart, acEnd - acStart));
1115: PetscCall(PetscViewerASCIISubtractTab(viewer, tablevel));
1116: PetscCall(PetscViewerPopFormat(viewer));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*@C
1121: DMAdaptorMonitorError - Prints the error norm at each iteration of an adaptation loop.
1123: Collective
1125: Input Parameters:
1126: + adaptor - the `DMAdaptor`
1127: . n - iteration number
1128: . odm - the original `DM`
1129: . adm - the adapted `DM`
1130: . Nf - number of fields
1131: . enorms - 2-norm error values for each field (may be estimated).
1132: . error - `Vec` of cellwise errors
1133: - vf - The viewer context
1135: Options Database Key:
1136: . -adaptor_monitor_error - Activates `DMAdaptorMonitorError()`
1138: Level: intermediate
1140: Note:
1141: This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
1142: to be used during the adaptation loop.
1144: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorErrorDrawLG()`
1145: @*/
1146: PetscErrorCode DMAdaptorMonitorError(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
1147: {
1148: PetscViewer viewer = vf->viewer;
1149: PetscViewerFormat format = vf->format;
1150: PetscInt tablevel, cStart, cEnd, acStart, acEnd;
1151: const char *prefix;
1153: PetscFunctionBegin;
1155: PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
1156: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
1158: PetscCall(PetscViewerPushFormat(viewer, format));
1159: PetscCall(PetscViewerASCIIAddTab(viewer, tablevel));
1160: if (n == 0 && prefix) PetscCall(PetscViewerASCIIPrintf(viewer, " Error norms for %s adaptation.\n", prefix));
1161: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " DMAdaptor Error norm %s", n, Nf > 1 ? "[" : ""));
1162: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1163: for (PetscInt f = 0; f < Nf; ++f) {
1164: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
1165: PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)enorms[f]));
1166: }
1167: PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
1168: PetscCall(DMPlexGetHeightStratum(adm, 0, &acStart, &acEnd));
1169: PetscCall(PetscViewerASCIIPrintf(viewer, " N: %" PetscInt_FMT " Nadapt: %" PetscInt_FMT "\n", cEnd - cStart, acEnd - acStart));
1170: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1171: PetscCall(PetscViewerASCIISubtractTab(viewer, tablevel));
1172: PetscCall(PetscViewerPopFormat(viewer));
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: /*@C
1177: DMAdaptorMonitorErrorDraw - Plots the error at each iteration of an iterative solver.
1179: Collective
1181: Input Parameters:
1182: + adaptor - the `DMAdaptor`
1183: . n - iteration number
1184: . odm - the original `DM`
1185: . adm - the adapted `DM`
1186: . Nf - number of fields
1187: . enorms - 2-norm error values for each field (may be estimated).
1188: . error - `Vec` of cellwise errors
1189: - vf - The viewer context
1191: Options Database Key:
1192: . -adaptor_monitor_error draw - Activates `DMAdaptorMonitorErrorDraw()`
1194: Level: intermediate
1196: Note:
1197: This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
1198: to be used during the adaptation loop.
1200: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDrawLG()`
1201: @*/
1202: PetscErrorCode DMAdaptorMonitorErrorDraw(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
1203: {
1204: PetscViewer viewer = vf->viewer;
1205: PetscViewerFormat format = vf->format;
1207: PetscFunctionBegin;
1209: PetscCall(PetscViewerPushFormat(viewer, format));
1210: PetscCall(PetscObjectSetName((PetscObject)error, "Error Estimator"));
1211: PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", (PetscObject)adaptor));
1212: PetscCall(VecView(error, viewer));
1213: PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", NULL));
1214: PetscCall(PetscViewerPopFormat(viewer));
1215: PetscFunctionReturn(PETSC_SUCCESS);
1216: }
1218: /*@C
1219: DMAdaptorMonitorErrorDrawLGCreate - Creates the context for the error plotter `DMAdaptorMonitorErrorDrawLG()`
1221: Collective
1223: Input Parameters:
1224: + viewer - The `PetscViewer`
1225: . format - The viewer format
1226: - ctx - An optional application context
1228: Output Parameter:
1229: . vf - The viewer context
1231: Level: intermediate
1233: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `PetscViewerMonitorGLSetUp()`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDrawLG()`
1234: @*/
1235: PetscErrorCode DMAdaptorMonitorErrorDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, PetscCtx ctx, PetscViewerAndFormat **vf)
1236: {
1237: DMAdaptor adaptor = (DMAdaptor)ctx;
1238: char **names;
1239: PetscInt Nf;
1241: PetscFunctionBegin;
1242: PetscCall(DMGetNumFields(adaptor->idm, &Nf));
1243: PetscCall(PetscMalloc1(Nf + 1, &names));
1244: for (PetscInt f = 0; f < Nf; ++f) {
1245: PetscObject disc;
1246: const char *fname;
1247: char lname[PETSC_MAX_PATH_LEN];
1249: PetscCall(DMGetField(adaptor->idm, f, NULL, &disc));
1250: PetscCall(PetscObjectGetName(disc, &fname));
1251: PetscCall(PetscStrncpy(lname, fname, PETSC_MAX_PATH_LEN));
1252: PetscCall(PetscStrlcat(lname, " Error", PETSC_MAX_PATH_LEN));
1253: PetscCall(PetscStrallocpy(lname, &names[f]));
1254: }
1255: PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
1256: (*vf)->data = ctx;
1257: PetscCall(PetscViewerMonitorLGSetUp(viewer, NULL, NULL, "Log Error Norm", Nf, (const char **)names, PETSC_DECIDE, PETSC_DECIDE, 400, 300));
1258: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(names[f]));
1259: PetscCall(PetscFree(names));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*@C
1264: DMAdaptorMonitorErrorDrawLG - Plots the error norm at each iteration of an adaptive loop.
1266: Collective
1268: Input Parameters:
1269: + adaptor - the `DMAdaptor`
1270: . n - iteration number
1271: . odm - the original `DM`
1272: . adm - the adapted `DM`
1273: . Nf - number of fields
1274: . enorms - 2-norm error values for each field (may be estimated).
1275: . error - `Vec` of cellwise errors
1276: - vf - The viewer context, obtained via `DMAdaptorMonitorErrorDrawLGCreate()`
1278: Options Database Key:
1279: . -adaptor_error draw::draw_lg - Activates `DMAdaptorMonitorErrorDrawLG()`
1281: Level: intermediate
1283: Notes:
1284: This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
1285: to be used during the adaptation loop.
1287: Call `DMAdaptorMonitorErrorDrawLGCreate()` to create the context needed for this monitor
1289: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorError()`,
1290: `DMAdaptorMonitorTrueResidualDrawLGCreate()`
1291: @*/
1292: PetscErrorCode DMAdaptorMonitorErrorDrawLG(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
1293: {
1294: PetscViewer viewer = vf->viewer;
1295: PetscViewerFormat format = vf->format;
1296: PetscDrawLG lg;
1297: PetscReal *x, *e;
1299: PetscFunctionBegin;
1301: PetscCall(PetscViewerDrawGetDrawLG(viewer, 0, &lg));
1302: PetscCall(PetscCalloc2(Nf, &x, Nf, &e));
1303: PetscCall(PetscViewerPushFormat(viewer, format));
1304: if (!n) PetscCall(PetscDrawLGReset(lg));
1305: for (PetscInt f = 0; f < Nf; ++f) {
1306: x[f] = (PetscReal)n;
1307: e[f] = enorms[f] > 0.0 ? PetscLog10Real(enorms[f]) : -15.;
1308: }
1309: PetscCall(PetscDrawLGAddPoint(lg, x, e));
1310: PetscCall(PetscDrawLGDraw(lg));
1311: PetscCall(PetscDrawLGSave(lg));
1312: PetscCall(PetscViewerPopFormat(viewer));
1313: PetscCall(PetscFree2(x, e));
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }
1317: /*@C
1318: DMAdaptorMonitorRegisterAll - Registers all of the mesh adaptation monitors in the `SNES` package.
1320: Not Collective
1322: Level: advanced
1324: .seealso: [](ch_snes), `SNES`, `DM`, `DMAdaptorMonitorRegister()`, `DMAdaptorRegister()`
1325: @*/
1326: PetscErrorCode DMAdaptorMonitorRegisterAll(void)
1327: {
1328: PetscFunctionBegin;
1329: if (DMAdaptorMonitorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1330: DMAdaptorMonitorRegisterAllCalled = PETSC_TRUE;
1332: PetscCall(DMAdaptorMonitorRegister("size", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorSize, NULL, NULL));
1333: PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorError, NULL, NULL));
1334: PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorErrorDraw, NULL, NULL));
1335: PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DRAW_LG, DMAdaptorMonitorErrorDrawLG, DMAdaptorMonitorErrorDrawLGCreate, NULL));
1336: PetscFunctionReturn(PETSC_SUCCESS);
1337: }
1339: static void identity(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[])
1340: {
1341: const PetscInt Nc = uOff[1] - uOff[0];
1343: for (PetscInt i = 0; i < Nc; ++i) f[i] = u[i];
1344: }
1346: 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[])
1347: {
1348: for (PetscInt i = 0; i < dim; ++i) {
1349: for (PetscInt j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
1350: }
1351: }
1353: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
1354: {
1355: PetscDS ds;
1356: PetscReal errorNorm = 0.;
1357: PetscInt numAdapt = adaptor->numSeq, adaptIter;
1358: PetscInt dim = 0, coordDim = 0, Nf = 0;
1359: void *ctx = NULL;
1360: MPI_Comm comm;
1362: PetscFunctionBegin;
1363: PetscCall(DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view"));
1364: PetscCall(VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view"));
1365: PetscCall(PetscObjectGetComm((PetscObject)adaptor, &comm));
1366: PetscCall(DMGetApplicationContext(adaptor->idm, &ctx));
1367: PetscCall(DMGetDS(adaptor->idm, &ds));
1368: PetscCall(PetscDSGetNumFields(ds, &Nf));
1369: if (adaptor->adaptCriterion != DM_ADAPTATION_REFINE) {
1370: PetscCall(DMGetDimension(adaptor->idm, &dim));
1371: PetscCall(DMGetCoordinateDim(adaptor->idm, &coordDim));
1372: PetscCheck(Nf != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot adapt with no fields present!");
1373: }
1375: /* Adapt until nothing changes */
1376: /* Adapt for a specified number of iterates */
1377: for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm)));
1378: for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
1379: PetscBool adapted = PETSC_FALSE;
1380: DM dm = adaptIter ? *adm : adaptor->idm, odm;
1381: Vec x = adaptIter ? *ax : inx, locX = NULL, ox;
1382: Vec error = NULL;
1384: if (adaptor->adaptCriterion != DM_ADAPTATION_REFINE) {
1385: PetscCall(DMGetLocalVector(dm, &locX));
1386: PetscCall(DMAdaptorPreAdapt(adaptor, locX));
1387: }
1388: if (doSolve) {
1389: SNES snes;
1391: PetscCall(DMAdaptorGetSolver(adaptor, &snes));
1392: PetscCall(SNESSolve(snes, NULL, x));
1393: }
1394: if (adaptor->adaptCriterion != DM_ADAPTATION_REFINE) {
1395: PetscCall(DMGlobalToLocalBegin(dm, x, INSERT_VALUES, locX));
1396: PetscCall(DMGlobalToLocalEnd(dm, x, INSERT_VALUES, locX));
1397: }
1398: PetscCall(VecViewFromOptions(x, (PetscObject)adaptor, "-adapt_primal_sol_vec_view"));
1399: switch (adaptor->adaptCriterion) {
1400: case DM_ADAPTATION_REFINE:
1401: PetscCall(DMRefine(dm, comm, &odm));
1402: PetscCheck(odm, comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
1403: adapted = PETSC_TRUE;
1404: PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, dm, 1, &errorNorm, NULL));
1405: break;
1406: case DM_ADAPTATION_LABEL: {
1407: /* Adapt DM
1408: Create local solution
1409: Reconstruct gradients (FVM) or solve adjoint equation (FEM)
1410: Produce cellwise error indicator */
1411: DM edm, plex;
1412: PetscDS ds;
1413: PetscFE efe;
1414: DMLabel adaptLabel;
1415: IS refineIS, coarsenIS;
1416: DMPolytopeType ct;
1417: PetscScalar errorVal;
1418: PetscInt nRefine, nCoarsen, cStart;
1420: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
1422: // TODO Move this creation to PreAdapt
1423: PetscCall(DMClone(dm, &edm));
1424: PetscCall(DMConvert(edm, DMPLEX, &plex));
1425: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, NULL));
1426: PetscCall(DMPlexGetCellType(plex, cStart, &ct));
1427: PetscCall(DMDestroy(&plex));
1428: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &efe));
1429: PetscCall(PetscObjectSetName((PetscObject)efe, "Error"));
1430: PetscCall(DMSetField(edm, 0, NULL, (PetscObject)efe));
1431: PetscCall(PetscFEDestroy(&efe));
1432: PetscCall(DMCreateDS(edm));
1433: PetscCall(DMGetGlobalVector(edm, &error));
1434: PetscCall(PetscObjectSetName((PetscObject)error, "Error Estimator"));
1436: PetscUseTypeMethod(adaptor, computeerrorindicator, locX, error);
1437: PetscCall(VecViewFromOptions(error, (PetscObject)adaptor, "-adapt_error_vec_view"));
1438: PetscCall(DMGetDS(edm, &ds));
1439: PetscCall(PetscDSSetObjective(ds, 0, identity));
1440: PetscCall(DMPlexComputeIntegralFEM(edm, error, &errorVal, NULL));
1441: errorNorm = PetscRealPart(errorVal);
1443: // Compute IS from VecTagger
1444: PetscCall(VecTaggerComputeIS(adaptor->refineTag, error, &refineIS, NULL));
1445: PetscCall(VecTaggerComputeIS(adaptor->coarsenTag, error, &coarsenIS, NULL));
1446: PetscCall(ISViewFromOptions(refineIS, (PetscObject)adaptor->refineTag, "-is_view"));
1447: PetscCall(ISViewFromOptions(coarsenIS, (PetscObject)adaptor->coarsenTag, "-is_view"));
1448: PetscCall(ISGetSize(refineIS, &nRefine));
1449: PetscCall(ISGetSize(coarsenIS, &nCoarsen));
1450: PetscCall(PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen));
1451: if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
1452: if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
1453: PetscCall(ISDestroy(&coarsenIS));
1454: PetscCall(ISDestroy(&refineIS));
1455: // Adapt DM from label
1456: if (nRefine || nCoarsen) {
1457: char oprefix[PETSC_MAX_PATH_LEN];
1458: const char *p;
1459: PetscBool flg;
1461: PetscCall(PetscOptionsHasName(NULL, adaptor->hdr.prefix, "-adapt_vec_view", &flg));
1462: if (flg) {
1463: Vec ref;
1465: PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
1466: PetscCall(VecViewFromOptions(ref, (PetscObject)adaptor, "-adapt_vec_view"));
1467: PetscCall(VecDestroy(&ref));
1468: }
1470: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &p));
1471: PetscCall(PetscStrncpy(oprefix, p, PETSC_MAX_PATH_LEN));
1472: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)dm, "adapt_"));
1473: PetscCall(DMAdaptLabel(dm, adaptLabel, &odm));
1474: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, oprefix));
1475: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)odm, oprefix));
1476: PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, odm, 1, &errorNorm, error));
1477: adapted = PETSC_TRUE;
1478: } else {
1479: PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, dm, 1, &errorNorm, error));
1480: }
1481: PetscCall(DMLabelDestroy(&adaptLabel));
1482: PetscCall(DMRestoreGlobalVector(edm, &error));
1483: PetscCall(DMDestroy(&edm));
1484: } break;
1485: case DM_ADAPTATION_METRIC: {
1486: DM dmGrad, dmHess, dmMetric, dmDet;
1487: Vec xGrad, xHess, metric, determinant;
1488: PetscReal N;
1489: DMLabel bdLabel = NULL, rgLabel = NULL;
1490: PetscBool higherOrder = PETSC_FALSE;
1491: PetscInt Nd = coordDim * coordDim, f, vStart, vEnd;
1492: 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[]);
1494: PetscCall(PetscMalloc(1, &funcs));
1495: funcs[0] = identityFunc;
1497: /* Setup finite element spaces */
1498: PetscCall(DMClone(dm, &dmGrad));
1499: PetscCall(DMClone(dm, &dmHess));
1500: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple fields not yet considered"); // TODO
1501: for (f = 0; f < Nf; ++f) {
1502: PetscFE fe, feGrad, feHess;
1503: PetscDualSpace Q;
1504: PetscSpace space;
1505: DM K;
1506: PetscQuadrature q;
1507: PetscInt Nc, qorder, p;
1508: const char *prefix;
1510: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
1511: PetscCall(PetscFEGetNumComponents(fe, &Nc));
1512: PetscCheck(Nc <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple components not yet considered"); // TODO
1513: PetscCall(PetscFEGetBasisSpace(fe, &space));
1514: PetscCall(PetscSpaceGetDegree(space, NULL, &p));
1515: if (p > 1) higherOrder = PETSC_TRUE;
1516: PetscCall(PetscFEGetDualSpace(fe, &Q));
1517: PetscCall(PetscDualSpaceGetDM(Q, &K));
1518: PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
1519: PetscCall(PetscFEGetQuadrature(fe, &q));
1520: PetscCall(PetscQuadratureGetOrder(q, &qorder));
1521: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix));
1522: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad));
1523: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess));
1524: PetscCall(DMSetField(dmGrad, f, NULL, (PetscObject)feGrad));
1525: PetscCall(DMSetField(dmHess, f, NULL, (PetscObject)feHess));
1526: PetscCall(DMCreateDS(dmGrad));
1527: PetscCall(DMCreateDS(dmHess));
1528: PetscCall(PetscFEDestroy(&feGrad));
1529: PetscCall(PetscFEDestroy(&feHess));
1530: }
1531: /* Compute vertexwise gradients from cellwise gradients */
1532: PetscCall(DMCreateLocalVector(dmGrad, &xGrad));
1533: PetscCall(VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view"));
1534: PetscCall(DMPlexComputeGradientClementInterpolant(dm, locX, xGrad));
1535: PetscCall(VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view"));
1536: /* Compute vertexwise Hessians from cellwise Hessians */
1537: PetscCall(DMCreateLocalVector(dmHess, &xHess));
1538: PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess));
1539: PetscCall(VecViewFromOptions(xHess, NULL, "-adapt_hessian_view"));
1540: PetscCall(VecDestroy(&xGrad));
1541: PetscCall(DMDestroy(&dmGrad));
1542: /* Compute L-p normalized metric */
1543: PetscCall(DMClone(dm, &dmMetric));
1544: N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
1545: // TODO This was where the old monitor was, figure out how to show metric and target N
1546: PetscCall(DMPlexMetricSetTargetComplexity(dmMetric, N));
1547: if (higherOrder) {
1548: /* Project Hessian into P1 space, if required */
1549: PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
1550: PetscCall(DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric));
1551: PetscCall(VecDestroy(&xHess));
1552: xHess = metric;
1553: }
1554: PetscCall(PetscFree(funcs));
1555: PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
1556: PetscCall(DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet));
1557: PetscCall(DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant));
1558: PetscCall(VecDestroy(&determinant));
1559: PetscCall(DMDestroy(&dmDet));
1560: PetscCall(VecDestroy(&xHess));
1561: PetscCall(DMDestroy(&dmHess));
1562: /* Adapt DM from metric */
1563: PetscCall(DMGetLabel(dm, "marker", &bdLabel));
1564: PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm));
1565: PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, odm, 1, &errorNorm, NULL));
1566: adapted = PETSC_TRUE;
1567: /* Cleanup */
1568: PetscCall(VecDestroy(&metric));
1569: PetscCall(DMDestroy(&dmMetric));
1570: } break;
1571: default:
1572: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
1573: }
1574: if (adaptor->adaptCriterion != DM_ADAPTATION_REFINE) {
1575: PetscCall(DMAdaptorPostAdapt(adaptor));
1576: PetscCall(DMRestoreLocalVector(dm, &locX));
1577: }
1578: /* If DM was adapted, replace objects and recreate solution */
1579: if (adapted) {
1580: const char *name;
1582: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1583: PetscCall(PetscObjectSetName((PetscObject)odm, name));
1584: /* Reconfigure solver */
1585: PetscCall(SNESReset(adaptor->snes));
1586: PetscCall(SNESSetDM(adaptor->snes, odm));
1587: PetscCall(DMAdaptorSetSolver(adaptor, adaptor->snes));
1588: if (Nf) PetscCall(DMPlexSetSNESLocalFEM(odm, PETSC_FALSE, ctx));
1589: PetscCall(SNESSetFromOptions(adaptor->snes));
1590: /* Transfer system */
1591: if (Nf) PetscCall(DMCopyDisc(dm, odm));
1592: /* Transfer solution */
1593: PetscCall(DMCreateGlobalVector(odm, &ox));
1594: PetscCall(PetscObjectGetName((PetscObject)x, &name));
1595: PetscCall(PetscObjectSetName((PetscObject)ox, name));
1596: PetscCall(DMAdaptorTransferSolution(adaptor, dm, x, odm, ox));
1597: /* Cleanup adaptivity info */
1598: if (adaptIter > 0) PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm)));
1599: PetscCall(DMForestSetAdaptivityForest(dm, NULL)); /* clear internal references to the previous dm */
1600: PetscCall(DMDestroy(&dm));
1601: PetscCall(VecDestroy(&x));
1602: *adm = odm;
1603: *ax = ox;
1604: } else {
1605: *adm = dm;
1606: *ax = x;
1607: adaptIter = numAdapt;
1608: }
1609: if (adaptIter < numAdapt - 1) {
1610: PetscCall(DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view"));
1611: PetscCall(VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view"));
1612: }
1613: }
1614: PetscCall(DMViewFromOptions(*adm, NULL, "-dm_adapt_view"));
1615: PetscCall(VecViewFromOptions(*ax, NULL, "-sol_adapt_view"));
1616: PetscFunctionReturn(PETSC_SUCCESS);
1617: }
1619: /*@
1620: DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem
1622: Not Collective
1624: Input Parameters:
1625: + adaptor - The `DMAdaptor` object
1626: . x - The global approximate solution
1627: - strategy - The adaptation strategy, see `DMAdaptationStrategy`
1629: Output Parameters:
1630: + adm - The adapted `DM`
1631: - ax - The adapted solution
1633: Options Database Keys:
1634: + -snes_adapt (initial|sequential|multigrid) - adaption strategy, see `DMAdaptationStrategy`
1635: . -adapt_gradient_view - View the Clement interpolant of the solution gradient
1636: . -adapt_hessian_view - View the Clement interpolant of the solution Hessian
1637: - -adapt_metric_view - View the metric tensor for adaptive mesh refinement
1639: Level: intermediate
1641: Note:
1642: When the mesh is adapted, one reference to `x` and one reference to the `DM` of the solver are
1643: consumed, matching the use in `SNESSolve()` grid sequencing where those objects are replaced by
1644: the adapted ones. A caller that keeps using the input objects must take an additional reference
1645: to each of them before calling this function.
1647: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptationStrategy`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`
1648: @*/
1649: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
1650: {
1651: PetscFunctionBegin;
1652: switch (strategy) {
1653: case DM_ADAPTATION_INITIAL:
1654: PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
1655: break;
1656: case DM_ADAPTATION_SEQUENTIAL:
1657: PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
1658: break;
1659: default:
1660: SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
1661: }
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*@C
1666: DMAdaptorGetMixedSetupFunction - Get the function setting up the mixed problem, if it exists
1668: Not Collective
1670: Input Parameter:
1671: . adaptor - the `DMAdaptor`
1673: Output Parameter:
1674: . setupFunc - the function setting up the mixed problem, or `NULL`
1676: Level: advanced
1678: .seealso: `DMAdaptor`, `DMAdaptorSetMixedSetupFunction()`, `DMAdaptorAdapt()`
1679: @*/
1680: PetscErrorCode DMAdaptorGetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (**setupFunc)(DMAdaptor, DM))
1681: {
1682: PetscFunctionBegin;
1684: PetscAssertPointer(setupFunc, 2);
1685: *setupFunc = adaptor->ops->mixedsetup;
1686: PetscFunctionReturn(PETSC_SUCCESS);
1687: }
1689: /*@C
1690: DMAdaptorSetMixedSetupFunction - Set the function setting up the mixed problem
1692: Not Collective
1694: Input Parameters:
1695: + adaptor - the `DMAdaptor`
1696: - setupFunc - the function setting up the mixed problem
1698: Calling sequence of setupFunc:
1699: + adaptor - the `DMAdaptor`
1700: - dm - the `DM`
1702: Level: advanced
1704: .seealso: `DMAdaptor`, `DMAdaptorGetMixedSetupFunction()`, `DMAdaptorAdapt()`
1705: @*/
1706: PetscErrorCode DMAdaptorSetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (*setupFunc)(DMAdaptor adaptor, DM dm))
1707: {
1708: PetscFunctionBegin;
1711: adaptor->ops->mixedsetup = setupFunc;
1712: PetscFunctionReturn(PETSC_SUCCESS);
1713: }
1715: /*@
1716: DMAdaptorGetCriterion - Get the adaptation criterion
1718: Not Collective
1720: Input Parameter:
1721: . adaptor - the `DMAdaptor`
1723: Output Parameter:
1724: . criterion - the criterion for adaptation
1726: Level: advanced
1728: .seealso: `DMAdaptor`, `DMAdaptorSetCriterion()`, `DMAdaptationCriterion`
1729: @*/
1730: PetscErrorCode DMAdaptorGetCriterion(DMAdaptor adaptor, DMAdaptationCriterion *criterion)
1731: {
1732: PetscFunctionBegin;
1734: PetscAssertPointer(criterion, 2);
1735: *criterion = adaptor->adaptCriterion;
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: /*@
1740: DMAdaptorSetCriterion - Set the adaptation criterion
1742: Not Collective
1744: Input Parameters:
1745: + adaptor - the `DMAdaptor`
1746: - criterion - the adaptation criterion
1748: Level: advanced
1750: .seealso: `DMAdaptor`, `DMAdaptorGetCriterion()`, `DMAdaptationCriterion`
1751: @*/
1752: PetscErrorCode DMAdaptorSetCriterion(DMAdaptor adaptor, DMAdaptationCriterion criterion)
1753: {
1754: PetscFunctionBegin;
1756: adaptor->adaptCriterion = criterion;
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: static PetscErrorCode DMAdaptorInitialize_Gradient(DMAdaptor adaptor)
1761: {
1762: PetscFunctionBegin;
1763: adaptor->ops->computeerrorindicator = DMAdaptorComputeErrorIndicator_Gradient;
1764: adaptor->ops->computecellerrorindicator = DMAdaptorComputeCellErrorIndicator_Gradient;
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor adaptor)
1769: {
1770: PetscFunctionBegin;
1772: adaptor->data = NULL;
1774: PetscCall(DMAdaptorInitialize_Gradient(adaptor));
1775: PetscFunctionReturn(PETSC_SUCCESS);
1776: }
1778: static PetscErrorCode DMAdaptorInitialize_Flux(DMAdaptor adaptor)
1779: {
1780: PetscFunctionBegin;
1781: adaptor->ops->computeerrorindicator = DMAdaptorComputeErrorIndicator_Flux;
1782: PetscFunctionReturn(PETSC_SUCCESS);
1783: }
1785: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor adaptor)
1786: {
1787: PetscFunctionBegin;
1789: adaptor->data = NULL;
1791: PetscCall(DMAdaptorInitialize_Flux(adaptor));
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: }