Actual source code: fv.c
1: #include <petsc/private/petscfvimpl.h>
2: #include <petscdmplex.h>
3: #include <petscdmplextransform.h>
4: #include <petscds.h>
6: PetscClassId PETSCLIMITER_CLASSID = 0;
8: PetscFunctionList PetscLimiterList = NULL;
9: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
11: PetscBool Limitercite = PETSC_FALSE;
12: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
13: " title = {Analysis of slope limiters on irregular grids},\n"
14: " journal = {AIAA paper},\n"
15: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
16: " volume = {490},\n"
17: " year = {2005}\n}\n";
19: /*@
20: PetscLimiterRegister - Adds a new `PetscLimiter` implementation
22: Not Collective, No Fortran Support
24: Input Parameters:
25: + sname - The name of a new user-defined creation routine
26: - function - The creation routine
28: Example Usage:
29: .vb
30: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
31: .ve
33: Then, your `PetscLimiter` type can be chosen with the procedural interface via
34: .vb
35: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
36: PetscLimiterSetType(PetscLimiter, "my_lim");
37: .ve
38: or at runtime via the option
39: .vb
40: -petsclimiter_type my_lim
41: .ve
43: Level: advanced
45: Note:
46: `PetscLimiterRegister()` may be called multiple times to add several user-defined PetscLimiters
48: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterRegisterAll()`
49: @*/
50: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
51: {
52: PetscFunctionBegin;
53: PetscCall(PetscFunctionListAdd(&PetscLimiterList, sname, function));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: /*@
58: PetscLimiterSetType - Builds a `PetscLimiter` for a given `PetscLimiterType`
60: Collective
62: Input Parameters:
63: + lim - The `PetscLimiter` object
64: - name - The kind of limiter
66: Options Database Key:
67: . -petsclimiter_type type - Sets the PetscLimiter type; use -help for a list of available types
69: Level: intermediate
71: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
72: @*/
73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
74: {
75: PetscErrorCode (*r)(PetscLimiter);
76: PetscBool match;
78: PetscFunctionBegin;
80: PetscCall(PetscObjectTypeCompare((PetscObject)lim, name, &match));
81: if (match) PetscFunctionReturn(PETSC_SUCCESS);
83: PetscCall(PetscLimiterRegisterAll());
84: PetscCall(PetscFunctionListFind(PetscLimiterList, name, &r));
85: PetscCheck(r, PetscObjectComm((PetscObject)lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
87: PetscTryTypeMethod(lim, destroy);
88: lim->ops->destroy = NULL;
90: PetscCall((*r)(lim));
91: PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@
96: PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.
98: Not Collective
100: Input Parameter:
101: . lim - The `PetscLimiter`
103: Output Parameter:
104: . name - The `PetscLimiterType`
106: Level: intermediate
108: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
112: PetscFunctionBegin;
114: PetscAssertPointer(name, 2);
115: PetscCall(PetscLimiterRegisterAll());
116: *name = ((PetscObject)lim)->type_name;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database
123: Collective
125: Input Parameters:
126: + A - the `PetscLimiter` object to view
127: . obj - optional object that provides the options prefix to use, pass `NULL` to use the options prefix of `A`
128: - name - command line option name
130: Options Database Key:
131: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
133: Level: intermediate
135: Note:
136: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
137: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
139: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`, `PetscOptionsCreateViewer()`
140: @*/
141: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
142: {
143: PetscFunctionBegin;
145: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: PetscLimiterView - Views a `PetscLimiter`
152: Collective
154: Input Parameters:
155: + lim - the `PetscLimiter` object to view
156: - v - the viewer
158: Level: beginner
160: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
161: @*/
162: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
163: {
164: PetscFunctionBegin;
166: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
167: PetscTryTypeMethod(lim, view, v);
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@
172: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
174: Collective
176: Input Parameter:
177: . lim - the `PetscLimiter` object to set options for
179: Level: intermediate
181: .seealso: `PetscLimiter`, `PetscLimiterView()`
182: @*/
183: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
184: {
185: const char *defaultType;
186: char name[256];
187: PetscBool flg;
189: PetscFunctionBegin;
191: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
192: else defaultType = ((PetscObject)lim)->type_name;
193: PetscCall(PetscLimiterRegisterAll());
195: PetscObjectOptionsBegin((PetscObject)lim);
196: PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, sizeof(name), &flg));
197: if (flg) {
198: PetscCall(PetscLimiterSetType(lim, name));
199: } else if (!((PetscObject)lim)->type_name) {
200: PetscCall(PetscLimiterSetType(lim, defaultType));
201: }
202: PetscTryTypeMethod(lim, setfromoptions);
203: /* process any options handlers added with PetscObjectAddOptionsHandler() */
204: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
205: PetscOptionsEnd();
206: PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
213: Collective
215: Input Parameter:
216: . lim - the `PetscLimiter` object to setup
218: Level: intermediate
220: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscLimiterDestroy()`
221: @*/
222: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
223: {
224: PetscFunctionBegin;
226: PetscTryTypeMethod(lim, setup);
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: PetscLimiterDestroy - Destroys a `PetscLimiter` object
233: Collective
235: Input Parameter:
236: . lim - the `PetscLimiter` object to destroy
238: Level: beginner
240: .seealso: `PetscLimiter`, `PetscLimiterView()`
241: @*/
242: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
243: {
244: PetscFunctionBegin;
245: if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);
248: if (--((PetscObject)*lim)->refct > 0) {
249: *lim = NULL;
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
252: ((PetscObject)*lim)->refct = 0;
254: PetscTryTypeMethod(*lim, destroy);
255: PetscCall(PetscHeaderDestroy(lim));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@
260: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
262: Collective
264: Input Parameter:
265: . comm - The communicator for the `PetscLimiter` object
267: Output Parameter:
268: . lim - The `PetscLimiter` object
270: Level: beginner
272: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
273: @*/
274: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
275: {
276: PetscLimiter l;
278: PetscFunctionBegin;
279: PetscAssertPointer(lim, 2);
280: PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
281: PetscCall(PetscFVInitializePackage());
283: PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));
285: *lim = l;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: PetscLimiterLimit - Limit the flux
292: Input Parameters:
293: + lim - The `PetscLimiter`
294: - flim - The input field
296: Output Parameter:
297: . phi - The limited field
299: Level: beginner
301: Note:
302: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
303: .vb
304: The classical flux-limited formulation is psi(r) where
306: r = (u[0] - u[-1]) / (u[1] - u[0])
308: The second order TVD region is bounded by
310: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
312: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
313: phi(r)(r+1)/2 in which the reconstructed interface values are
315: u(v) = u[0] + phi(r) (grad u)[0] v
317: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
319: phi_minmod(r) = 2 min(1/(1+r),r/(1+r)) phi_superbee(r) = 2 min(2/(1+r), 2r/(1+r), max(1,r)/(1+r))
321: For a nicer symmetric formulation, rewrite in terms of
323: f = (u[0] - u[-1]) / (u[1] - u[-1])
325: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
327: phi(r) = phi(1/r)
329: becomes
331: w(f) = w(1-f).
333: The limiters below implement this final form w(f). The reference methods are
335: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
336: .ve
338: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
339: @*/
340: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
341: {
342: PetscFunctionBegin;
344: PetscAssertPointer(phi, 3);
345: PetscUseTypeMethod(lim, limit, flim, phi);
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
350: {
351: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
353: PetscFunctionBegin;
354: PetscCall(PetscFree(l));
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
359: {
360: PetscViewerFormat format;
362: PetscFunctionBegin;
363: PetscCall(PetscViewerGetFormat(viewer, &format));
364: PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
369: {
370: PetscBool isascii;
372: PetscFunctionBegin;
375: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
376: if (isascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
381: {
382: PetscFunctionBegin;
383: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
388: {
389: PetscFunctionBegin;
390: lim->ops->view = PetscLimiterView_Sin;
391: lim->ops->destroy = PetscLimiterDestroy_Sin;
392: lim->ops->limit = PetscLimiterLimit_Sin;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*MC
397: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
399: Level: intermediate
401: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
402: M*/
404: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
405: {
406: PetscLimiter_Sin *l;
408: PetscFunctionBegin;
410: PetscCall(PetscNew(&l));
411: lim->data = l;
413: PetscCall(PetscLimiterInitialize_Sin(lim));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
418: {
419: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
421: PetscFunctionBegin;
422: PetscCall(PetscFree(l));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
427: {
428: PetscViewerFormat format;
430: PetscFunctionBegin;
431: PetscCall(PetscViewerGetFormat(viewer, &format));
432: PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
437: {
438: PetscBool isascii;
440: PetscFunctionBegin;
443: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
444: if (isascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
449: {
450: PetscFunctionBegin;
451: *phi = 0.0;
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
456: {
457: PetscFunctionBegin;
458: lim->ops->view = PetscLimiterView_Zero;
459: lim->ops->destroy = PetscLimiterDestroy_Zero;
460: lim->ops->limit = PetscLimiterLimit_Zero;
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*MC
465: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
467: Level: intermediate
469: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
470: M*/
472: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
473: {
474: PetscLimiter_Zero *l;
476: PetscFunctionBegin;
478: PetscCall(PetscNew(&l));
479: lim->data = l;
481: PetscCall(PetscLimiterInitialize_Zero(lim));
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
486: {
487: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
489: PetscFunctionBegin;
490: PetscCall(PetscFree(l));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
495: {
496: PetscViewerFormat format;
498: PetscFunctionBegin;
499: PetscCall(PetscViewerGetFormat(viewer, &format));
500: PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
505: {
506: PetscBool isascii;
508: PetscFunctionBegin;
511: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
512: if (isascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
517: {
518: PetscFunctionBegin;
519: *phi = 1.0;
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
524: {
525: PetscFunctionBegin;
526: lim->ops->view = PetscLimiterView_None;
527: lim->ops->destroy = PetscLimiterDestroy_None;
528: lim->ops->limit = PetscLimiterLimit_None;
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: /*MC
533: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
535: Level: intermediate
537: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
538: M*/
540: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
541: {
542: PetscLimiter_None *l;
544: PetscFunctionBegin;
546: PetscCall(PetscNew(&l));
547: lim->data = l;
549: PetscCall(PetscLimiterInitialize_None(lim));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
554: {
555: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
557: PetscFunctionBegin;
558: PetscCall(PetscFree(l));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
563: {
564: PetscViewerFormat format;
566: PetscFunctionBegin;
567: PetscCall(PetscViewerGetFormat(viewer, &format));
568: PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
573: {
574: PetscBool isascii;
576: PetscFunctionBegin;
579: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
580: if (isascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
585: {
586: PetscFunctionBegin;
587: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
592: {
593: PetscFunctionBegin;
594: lim->ops->view = PetscLimiterView_Minmod;
595: lim->ops->destroy = PetscLimiterDestroy_Minmod;
596: lim->ops->limit = PetscLimiterLimit_Minmod;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*MC
601: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
603: Level: intermediate
605: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
606: M*/
608: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
609: {
610: PetscLimiter_Minmod *l;
612: PetscFunctionBegin;
614: PetscCall(PetscNew(&l));
615: lim->data = l;
617: PetscCall(PetscLimiterInitialize_Minmod(lim));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
622: {
623: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
625: PetscFunctionBegin;
626: PetscCall(PetscFree(l));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
631: {
632: PetscViewerFormat format;
634: PetscFunctionBegin;
635: PetscCall(PetscViewerGetFormat(viewer, &format));
636: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
641: {
642: PetscBool isascii;
644: PetscFunctionBegin;
647: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
648: if (isascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
653: {
654: PetscFunctionBegin;
655: *phi = PetscMax(0, 4 * f * (1 - f));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
660: {
661: PetscFunctionBegin;
662: lim->ops->view = PetscLimiterView_VanLeer;
663: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
664: lim->ops->limit = PetscLimiterLimit_VanLeer;
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /*MC
669: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
671: Level: intermediate
673: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
674: M*/
676: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
677: {
678: PetscLimiter_VanLeer *l;
680: PetscFunctionBegin;
682: PetscCall(PetscNew(&l));
683: lim->data = l;
685: PetscCall(PetscLimiterInitialize_VanLeer(lim));
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
690: {
691: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
693: PetscFunctionBegin;
694: PetscCall(PetscFree(l));
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
699: {
700: PetscViewerFormat format;
702: PetscFunctionBegin;
703: PetscCall(PetscViewerGetFormat(viewer, &format));
704: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
709: {
710: PetscBool isascii;
712: PetscFunctionBegin;
715: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
716: if (isascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
721: {
722: PetscFunctionBegin;
723: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
728: {
729: PetscFunctionBegin;
730: lim->ops->view = PetscLimiterView_VanAlbada;
731: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
732: lim->ops->limit = PetscLimiterLimit_VanAlbada;
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*MC
737: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
739: Level: intermediate
741: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
742: M*/
744: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
745: {
746: PetscLimiter_VanAlbada *l;
748: PetscFunctionBegin;
750: PetscCall(PetscNew(&l));
751: lim->data = l;
753: PetscCall(PetscLimiterInitialize_VanAlbada(lim));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
758: {
759: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
761: PetscFunctionBegin;
762: PetscCall(PetscFree(l));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
767: {
768: PetscViewerFormat format;
770: PetscFunctionBegin;
771: PetscCall(PetscViewerGetFormat(viewer, &format));
772: PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
777: {
778: PetscBool isascii;
780: PetscFunctionBegin;
783: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
784: if (isascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
789: {
790: PetscFunctionBegin;
791: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }
795: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
796: {
797: PetscFunctionBegin;
798: lim->ops->view = PetscLimiterView_Superbee;
799: lim->ops->destroy = PetscLimiterDestroy_Superbee;
800: lim->ops->limit = PetscLimiterLimit_Superbee;
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*MC
805: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
807: Level: intermediate
809: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
810: M*/
812: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
813: {
814: PetscLimiter_Superbee *l;
816: PetscFunctionBegin;
818: PetscCall(PetscNew(&l));
819: lim->data = l;
821: PetscCall(PetscLimiterInitialize_Superbee(lim));
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
826: {
827: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
829: PetscFunctionBegin;
830: PetscCall(PetscFree(l));
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
835: {
836: PetscViewerFormat format;
838: PetscFunctionBegin;
839: PetscCall(PetscViewerGetFormat(viewer, &format));
840: PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
845: {
846: PetscBool isascii;
848: PetscFunctionBegin;
851: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
852: if (isascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: /* aka Barth-Jespersen */
857: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
858: {
859: PetscFunctionBegin;
860: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
865: {
866: PetscFunctionBegin;
867: lim->ops->view = PetscLimiterView_MC;
868: lim->ops->destroy = PetscLimiterDestroy_MC;
869: lim->ops->limit = PetscLimiterLimit_MC;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*MC
874: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
876: Level: intermediate
878: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
879: M*/
881: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
882: {
883: PetscLimiter_MC *l;
885: PetscFunctionBegin;
887: PetscCall(PetscNew(&l));
888: lim->data = l;
890: PetscCall(PetscLimiterInitialize_MC(lim));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: PetscClassId PETSCFV_CLASSID = 0;
896: PetscFunctionList PetscFVList = NULL;
897: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
899: /*@
900: PetscFVRegister - Adds a new `PetscFV` implementation
902: Not Collective, No Fortran Support
904: Input Parameters:
905: + sname - The name of a new user-defined creation routine
906: - function - The creation routine itself
908: Example Usage:
909: .vb
910: PetscFVRegister("my_fv", MyPetscFVCreate);
911: .ve
913: Then, your PetscFV type can be chosen with the procedural interface via
914: .vb
915: PetscFVCreate(MPI_Comm, PetscFV *);
916: PetscFVSetType(PetscFV, "my_fv");
917: .ve
918: or at runtime via the option
919: .vb
920: -petscfv_type my_fv
921: .ve
923: Level: advanced
925: Note:
926: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
928: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`
929: @*/
930: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
931: {
932: PetscFunctionBegin;
933: PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: /*@
938: PetscFVSetType - Builds a particular `PetscFV`
940: Collective
942: Input Parameters:
943: + fvm - The `PetscFV` object
944: - name - The type of FVM space
946: Options Database Key:
947: . -petscfv_type type - Sets the `PetscFVType`; use -help for a list of available types
949: Level: intermediate
951: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
952: @*/
953: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
954: {
955: PetscErrorCode (*r)(PetscFV);
956: PetscBool match;
958: PetscFunctionBegin;
960: PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
961: if (match) PetscFunctionReturn(PETSC_SUCCESS);
963: PetscCall(PetscFVRegisterAll());
964: PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
965: PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
967: PetscTryTypeMethod(fvm, destroy);
968: fvm->ops->destroy = NULL;
970: PetscCall((*r)(fvm));
971: PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
978: Not Collective
980: Input Parameter:
981: . fvm - The `PetscFV`
983: Output Parameter:
984: . name - The `PetscFVType` name
986: Level: intermediate
988: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
989: @*/
990: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
991: {
992: PetscFunctionBegin;
994: PetscAssertPointer(name, 2);
995: PetscCall(PetscFVRegisterAll());
996: *name = ((PetscObject)fvm)->type_name;
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*@
1001: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
1003: Collective
1005: Input Parameters:
1006: + A - the `PetscFV` object
1007: . obj - optional object that provides the options prefix, pass `NULL` to use the options prefix of `A`
1008: - name - command line option name
1010: Options Database Key:
1011: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
1013: Level: intermediate
1015: Note:
1016: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
1017: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
1019: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`, `PetscOptionsCreateViewer()`
1020: @*/
1021: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1022: {
1023: PetscFunctionBegin;
1025: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: /*@
1030: PetscFVView - Views a `PetscFV`
1032: Collective
1034: Input Parameters:
1035: + fvm - the `PetscFV` object to view
1036: - v - the viewer
1038: Level: beginner
1040: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1041: @*/
1042: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1043: {
1044: PetscFunctionBegin;
1046: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1047: PetscTryTypeMethod(fvm, view, v);
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
1054: Collective
1056: Input Parameter:
1057: . fvm - the `PetscFV` object to set options for
1059: Options Database Key:
1060: . -petscfv_compute_gradients (true|false) - Determines whether cell gradients are calculated
1062: Level: intermediate
1064: .seealso: `PetscFV`, `PetscFVView()`
1065: @*/
1066: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1067: {
1068: const char *defaultType;
1069: char name[256];
1070: PetscBool flg;
1072: PetscFunctionBegin;
1074: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1075: else defaultType = ((PetscObject)fvm)->type_name;
1076: PetscCall(PetscFVRegisterAll());
1078: PetscObjectOptionsBegin((PetscObject)fvm);
1079: PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, sizeof(name), &flg));
1080: if (flg) PetscCall(PetscFVSetType(fvm, name));
1081: else if (!((PetscObject)fvm)->type_name) PetscCall(PetscFVSetType(fvm, defaultType));
1082: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1083: PetscTryTypeMethod(fvm, setfromoptions);
1084: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1085: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1086: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1087: PetscOptionsEnd();
1088: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1089: PetscFunctionReturn(PETSC_SUCCESS);
1090: }
1092: /*@
1093: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1095: Collective
1097: Input Parameter:
1098: . fvm - the `PetscFV` object to setup
1100: Level: intermediate
1102: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1103: @*/
1104: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1105: {
1106: PetscFunctionBegin;
1108: PetscCall(PetscLimiterSetUp(fvm->limiter));
1109: PetscTryTypeMethod(fvm, setup);
1110: PetscFunctionReturn(PETSC_SUCCESS);
1111: }
1113: /*@
1114: PetscFVDestroy - Destroys a `PetscFV` object
1116: Collective
1118: Input Parameter:
1119: . fvm - the `PetscFV` object to destroy
1121: Level: beginner
1123: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1124: @*/
1125: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1126: {
1127: PetscFunctionBegin;
1128: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1131: if (--((PetscObject)*fvm)->refct > 0) {
1132: *fvm = NULL;
1133: PetscFunctionReturn(PETSC_SUCCESS);
1134: }
1135: ((PetscObject)*fvm)->refct = 0;
1137: for (PetscInt i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1138: PetscCall(PetscFree((*fvm)->componentNames));
1139: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1140: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1141: PetscCall(PetscFree((*fvm)->fluxWork));
1142: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1143: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1145: PetscTryTypeMethod(*fvm, destroy);
1146: PetscCall(PetscHeaderDestroy(fvm));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: /*@
1151: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1153: Collective
1155: Input Parameter:
1156: . comm - The communicator for the `PetscFV` object
1158: Output Parameter:
1159: . fvm - The `PetscFV` object
1161: Level: beginner
1163: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1164: @*/
1165: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1166: {
1167: PetscFV f;
1169: PetscFunctionBegin;
1170: PetscAssertPointer(fvm, 2);
1171: PetscCall(PetscFVInitializePackage());
1173: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1174: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1175: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1176: f->numComponents = 1;
1177: f->dim = 0;
1178: f->computeGradients = PETSC_FALSE;
1179: f->fluxWork = NULL;
1180: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1182: *fvm = f;
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: /*@
1187: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1189: Logically Collective
1191: Input Parameters:
1192: + fvm - the `PetscFV` object
1193: - lim - The `PetscLimiter`
1195: Level: intermediate
1197: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1198: @*/
1199: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1200: {
1201: PetscFunctionBegin;
1204: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1205: PetscCall(PetscObjectReference((PetscObject)lim));
1206: fvm->limiter = lim;
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@
1211: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1213: Not Collective
1215: Input Parameter:
1216: . fvm - the `PetscFV` object
1218: Output Parameter:
1219: . lim - The `PetscLimiter`
1221: Level: intermediate
1223: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1224: @*/
1225: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1226: {
1227: PetscFunctionBegin;
1229: PetscAssertPointer(lim, 2);
1230: *lim = fvm->limiter;
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: /*@
1235: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1237: Logically Collective
1239: Input Parameters:
1240: + fvm - the `PetscFV` object
1241: - comp - The number of components
1243: Level: intermediate
1245: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1246: @*/
1247: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1248: {
1249: PetscFunctionBegin;
1251: if (fvm->numComponents != comp) {
1252: for (PetscInt i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1253: PetscCall(PetscFree(fvm->componentNames));
1254: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1255: }
1256: fvm->numComponents = comp;
1257: PetscCall(PetscFree(fvm->fluxWork));
1258: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1259: PetscFunctionReturn(PETSC_SUCCESS);
1260: }
1262: /*@
1263: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1265: Not Collective
1267: Input Parameter:
1268: . fvm - the `PetscFV` object
1270: Output Parameter:
1271: . comp - The number of components
1273: Level: intermediate
1275: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1276: @*/
1277: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1278: {
1279: PetscFunctionBegin;
1281: PetscAssertPointer(comp, 2);
1282: *comp = fvm->numComponents;
1283: PetscFunctionReturn(PETSC_SUCCESS);
1284: }
1286: /*@
1287: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1289: Logically Collective
1291: Input Parameters:
1292: + fvm - the `PetscFV` object
1293: . comp - the component number
1294: - name - the component name
1296: Level: intermediate
1298: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1299: @*/
1300: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1301: {
1302: PetscFunctionBegin;
1303: PetscCall(PetscFree(fvm->componentNames[comp]));
1304: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /*@
1309: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1311: Logically Collective
1313: Input Parameters:
1314: + fvm - the `PetscFV` object
1315: - comp - the component number
1317: Output Parameter:
1318: . name - the component name
1320: Level: intermediate
1322: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1323: @*/
1324: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char *name[])
1325: {
1326: PetscFunctionBegin;
1327: *name = fvm->componentNames[comp];
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: /*@
1332: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1334: Logically Collective
1336: Input Parameters:
1337: + fvm - the `PetscFV` object
1338: - dim - The spatial dimension
1340: Level: intermediate
1342: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
1343: @*/
1344: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1345: {
1346: PetscFunctionBegin;
1348: fvm->dim = dim;
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /*@
1353: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1355: Not Collective
1357: Input Parameter:
1358: . fvm - the `PetscFV` object
1360: Output Parameter:
1361: . dim - The spatial dimension
1363: Level: intermediate
1365: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1366: @*/
1367: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1368: {
1369: PetscFunctionBegin;
1371: PetscAssertPointer(dim, 2);
1372: *dim = fvm->dim;
1373: PetscFunctionReturn(PETSC_SUCCESS);
1374: }
1376: /*@
1377: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1379: Logically Collective
1381: Input Parameters:
1382: + fvm - the `PetscFV` object
1383: - computeGradients - Flag to compute cell gradients
1385: Level: intermediate
1387: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1388: @*/
1389: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1390: {
1391: PetscFunctionBegin;
1393: fvm->computeGradients = computeGradients;
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: /*@
1398: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1400: Not Collective
1402: Input Parameter:
1403: . fvm - the `PetscFV` object
1405: Output Parameter:
1406: . computeGradients - Flag to compute cell gradients
1408: Level: intermediate
1410: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1411: @*/
1412: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1413: {
1414: PetscFunctionBegin;
1416: PetscAssertPointer(computeGradients, 2);
1417: *computeGradients = fvm->computeGradients;
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: /*@
1422: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1424: Logically Collective
1426: Input Parameters:
1427: + fvm - the `PetscFV` object
1428: - q - The `PetscQuadrature`
1430: Level: intermediate
1432: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1433: @*/
1434: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1435: {
1436: PetscFunctionBegin;
1438: PetscCall(PetscObjectReference((PetscObject)q));
1439: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1440: fvm->quadrature = q;
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /*@
1445: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1447: Not Collective
1449: Input Parameter:
1450: . fvm - the `PetscFV` object
1452: Output Parameter:
1453: . q - The `PetscQuadrature`
1455: Level: intermediate
1457: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1458: @*/
1459: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1460: {
1461: PetscFunctionBegin;
1463: PetscAssertPointer(q, 2);
1464: if (!fvm->quadrature) {
1465: /* Create default 1-point quadrature */
1466: PetscReal *points, *weights;
1468: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1469: PetscCall(PetscCalloc1(fvm->dim, &points));
1470: PetscCall(PetscMalloc1(1, &weights));
1471: weights[0] = 1.0;
1472: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1473: }
1474: *q = fvm->quadrature;
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: /*@
1479: PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`
1481: Not Collective
1483: Input Parameters:
1484: + fvm - The `PetscFV` object
1485: - ct - The `DMPolytopeType` for the cell
1487: Level: intermediate
1489: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1490: @*/
1491: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1492: {
1493: DM K;
1494: PetscInt dim, Nc;
1496: PetscFunctionBegin;
1497: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1498: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1499: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1500: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1501: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1502: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1503: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1504: PetscCall(DMDestroy(&K));
1505: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1506: // Should we be using PetscFVGetQuadrature() here?
1507: for (PetscInt c = 0; c < Nc; ++c) {
1508: PetscQuadrature qc;
1509: PetscReal *points, *weights;
1511: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1512: PetscCall(PetscCalloc1(dim, &points));
1513: PetscCall(PetscCalloc1(Nc, &weights));
1514: weights[c] = 1.0;
1515: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1516: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1517: PetscCall(PetscQuadratureDestroy(&qc));
1518: }
1519: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }
1523: /*@
1524: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1526: Not Collective
1528: Input Parameter:
1529: . fvm - The `PetscFV` object
1531: Output Parameter:
1532: . sp - The `PetscDualSpace` object
1534: Level: intermediate
1536: Developer Notes:
1537: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1539: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1540: @*/
1541: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1542: {
1543: PetscFunctionBegin;
1545: PetscAssertPointer(sp, 2);
1546: if (!fvm->dualSpace) {
1547: PetscInt dim;
1549: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1550: PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1551: }
1552: *sp = fvm->dualSpace;
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: /*@
1557: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1559: Not Collective
1561: Input Parameters:
1562: + fvm - The `PetscFV` object
1563: - sp - The `PetscDualSpace` object
1565: Level: intermediate
1567: Note:
1568: A simple dual space is provided automatically, and the user typically will not need to override it.
1570: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1571: @*/
1572: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1573: {
1574: PetscFunctionBegin;
1577: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1578: fvm->dualSpace = sp;
1579: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1586: Not Collective
1588: Input Parameter:
1589: . fvm - The `PetscFV` object
1591: Output Parameter:
1592: . T - The basis function values and derivatives at quadrature points
1594: Level: intermediate
1596: Note:
1597: .vb
1598: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1599: T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1600: T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1601: .ve
1603: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1604: @*/
1605: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1606: {
1607: PetscInt npoints;
1608: const PetscReal *points;
1610: PetscFunctionBegin;
1612: PetscAssertPointer(T, 2);
1613: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1614: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1615: *T = fvm->T;
1616: PetscFunctionReturn(PETSC_SUCCESS);
1617: }
1619: /*@
1620: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1622: Not Collective
1624: Input Parameters:
1625: + fvm - The `PetscFV` object
1626: . nrepl - The number of replicas
1627: . npoints - The number of tabulation points in a replica
1628: . points - The tabulation point coordinates
1629: - K - The order of derivative to tabulate
1631: Output Parameter:
1632: . T - The basis function values and derivative at tabulation points
1634: Level: intermediate
1636: Note:
1637: .vb
1638: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1639: T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1640: T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1641: .ve
1643: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1644: @*/
1645: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1646: {
1647: PetscInt pdim; // Dimension of approximation space P
1648: PetscInt cdim; // Spatial dimension
1649: PetscInt Nc; // Field components
1650: PetscInt k, p, d, c, e;
1652: PetscFunctionBegin;
1653: if (!npoints || K < 0) {
1654: *T = NULL;
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: PetscAssertPointer(points, 4);
1659: PetscAssertPointer(T, 6);
1660: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1661: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1662: pdim = Nc;
1663: PetscCall(PetscMalloc1(1, T));
1664: (*T)->K = !cdim ? 0 : K;
1665: (*T)->Nr = nrepl;
1666: (*T)->Np = npoints;
1667: (*T)->Nb = pdim;
1668: (*T)->Nc = Nc;
1669: (*T)->cdim = cdim;
1670: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1671: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1672: if (K >= 0) {
1673: for (p = 0; p < nrepl * npoints; ++p)
1674: for (d = 0; d < pdim; ++d)
1675: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1676: }
1677: if (K >= 1) {
1678: for (p = 0; p < nrepl * npoints; ++p)
1679: for (d = 0; d < pdim; ++d)
1680: for (c = 0; c < Nc; ++c)
1681: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1682: }
1683: if (K >= 2) {
1684: for (p = 0; p < nrepl * npoints; ++p)
1685: for (d = 0; d < pdim; ++d)
1686: for (c = 0; c < Nc; ++c)
1687: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1688: }
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1695: Input Parameters:
1696: + fvm - The `PetscFV` object
1697: . numFaces - The number of cell faces which are not constrained
1698: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1700: Output Parameter:
1701: . grad - the gradient
1703: Level: advanced
1705: .seealso: `PetscFV`, `PetscFVCreate()`
1706: @*/
1707: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1708: {
1709: PetscFunctionBegin;
1711: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1712: PetscFunctionReturn(PETSC_SUCCESS);
1713: }
1715: /*@
1716: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1718: Not Collective
1720: Input Parameters:
1721: + fvm - The `PetscFV` object for the field being integrated
1722: . prob - The `PetscDS` specifying the discretizations and continuum functions
1723: . field - The field being integrated
1724: . Nf - The number of faces in the chunk
1725: . fgeom - The face geometry for each face in the chunk
1726: . neighborVol - The volume for each pair of cells in the chunk
1727: . uL - The state from the cell on the left
1728: - uR - The state from the cell on the right
1730: Output Parameters:
1731: + fluxL - the left fluxes for each face
1732: - fluxR - the right fluxes for each face
1734: Level: developer
1736: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1737: @*/
1738: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1739: {
1740: PetscFunctionBegin;
1742: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: /*@
1747: PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.
1749: Input Parameter:
1750: . fv - The initial `PetscFV`
1752: Output Parameter:
1753: . fvNew - A clone of the `PetscFV`
1755: Level: advanced
1757: Notes:
1758: This is typically used to change the number of components.
1760: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1761: @*/
1762: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1763: {
1764: PetscDualSpace Q;
1765: DM K;
1766: PetscQuadrature q;
1767: PetscInt Nc, cdim;
1769: PetscFunctionBegin;
1770: PetscCall(PetscFVGetDualSpace(fv, &Q));
1771: PetscCall(PetscFVGetQuadrature(fv, &q));
1772: PetscCall(PetscDualSpaceGetDM(Q, &K));
1774: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1775: PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1776: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1777: PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1778: PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1779: PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1780: PetscCall(PetscFVSetQuadrature(*fvNew, q));
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: /*@
1785: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1786: smaller copies.
1788: Input Parameter:
1789: . fv - The initial `PetscFV`
1791: Output Parameter:
1792: . fvRef - The refined `PetscFV`
1794: Level: advanced
1796: Notes:
1797: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1798: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1799: interpolation between regularly refined meshes.
1801: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1802: @*/
1803: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1804: {
1805: PetscDualSpace Q, Qref;
1806: DM K, Kref;
1807: PetscQuadrature q, qref;
1808: DMPolytopeType ct;
1809: DMPlexTransform tr;
1810: PetscReal *v0;
1811: PetscReal *jac, *invjac;
1812: PetscInt numComp, numSubelements, s;
1814: PetscFunctionBegin;
1815: PetscCall(PetscFVGetDualSpace(fv, &Q));
1816: PetscCall(PetscFVGetQuadrature(fv, &q));
1817: PetscCall(PetscDualSpaceGetDM(Q, &K));
1818: /* Create dual space */
1819: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1820: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1821: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1822: PetscCall(DMDestroy(&Kref));
1823: PetscCall(PetscDualSpaceSetUp(Qref));
1824: /* Create volume */
1825: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1826: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1827: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1828: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1829: PetscCall(PetscFVSetUp(*fvRef));
1830: /* Create quadrature */
1831: PetscCall(DMPlexGetCellType(K, 0, &ct));
1832: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1833: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1834: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1835: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1836: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1837: for (s = 0; s < numSubelements; ++s) {
1838: PetscQuadrature qs;
1839: const PetscReal *points, *weights;
1840: PetscReal *p, *w;
1841: PetscInt dim, Nc, npoints, np;
1843: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1844: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1845: np = npoints / numSubelements;
1846: PetscCall(PetscMalloc1(np * dim, &p));
1847: PetscCall(PetscMalloc1(np * Nc, &w));
1848: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1849: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1850: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1851: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1852: PetscCall(PetscQuadratureDestroy(&qs));
1853: }
1854: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1855: PetscCall(DMPlexTransformDestroy(&tr));
1856: PetscCall(PetscQuadratureDestroy(&qref));
1857: PetscCall(PetscDualSpaceDestroy(&Qref));
1858: PetscFunctionReturn(PETSC_SUCCESS);
1859: }
1861: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1862: {
1863: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1865: PetscFunctionBegin;
1866: PetscCall(PetscFree(b));
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1871: {
1872: PetscInt Nc;
1873: PetscViewerFormat format;
1875: PetscFunctionBegin;
1876: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1877: PetscCall(PetscViewerGetFormat(viewer, &format));
1878: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1879: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1880: for (PetscInt c = 0; c < Nc; c++) {
1881: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1882: }
1883: PetscFunctionReturn(PETSC_SUCCESS);
1884: }
1886: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1887: {
1888: PetscBool isascii;
1890: PetscFunctionBegin;
1893: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1894: if (isascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1895: PetscFunctionReturn(PETSC_SUCCESS);
1896: }
1898: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1899: {
1900: PetscInt dim;
1902: PetscFunctionBegin;
1903: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1904: for (PetscInt f = 0; f < numFaces; ++f) {
1905: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1906: }
1907: PetscFunctionReturn(PETSC_SUCCESS);
1908: }
1910: /*
1911: neighborVol[f*2+0] contains the left geom
1912: neighborVol[f*2+1] contains the right geom
1913: */
1914: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1915: {
1916: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1917: void *rctx;
1918: PetscScalar *flux = fvm->fluxWork;
1919: const PetscScalar *constants;
1920: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1922: PetscFunctionBegin;
1923: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1924: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1925: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1926: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1927: PetscCall(PetscDSGetContext(prob, field, &rctx));
1928: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1929: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1930: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1931: for (f = 0; f < Nf; ++f) {
1932: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1933: for (d = 0; d < pdim; ++d) {
1934: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1935: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1936: }
1937: }
1938: PetscFunctionReturn(PETSC_SUCCESS);
1939: }
1941: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1942: {
1943: PetscFunctionBegin;
1944: fvm->ops->setfromoptions = NULL;
1945: fvm->ops->view = PetscFVView_Upwind;
1946: fvm->ops->destroy = PetscFVDestroy_Upwind;
1947: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1948: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: /*MC
1953: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1955: Level: intermediate
1957: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1958: M*/
1960: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1961: {
1962: PetscFV_Upwind *b;
1964: PetscFunctionBegin;
1966: PetscCall(PetscNew(&b));
1967: fvm->data = b;
1969: PetscCall(PetscFVInitialize_Upwind(fvm));
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: #include <petscblaslapack.h>
1975: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1976: {
1977: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1979: PetscFunctionBegin;
1980: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1981: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1982: PetscCall(PetscFree(ls));
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1986: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1987: {
1988: PetscInt Nc;
1989: PetscViewerFormat format;
1991: PetscFunctionBegin;
1992: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1993: PetscCall(PetscViewerGetFormat(viewer, &format));
1994: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1995: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1996: for (PetscInt c = 0; c < Nc; c++) {
1997: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1998: }
1999: PetscFunctionReturn(PETSC_SUCCESS);
2000: }
2002: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2003: {
2004: PetscBool isascii;
2006: PetscFunctionBegin;
2009: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2010: if (isascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2011: PetscFunctionReturn(PETSC_SUCCESS);
2012: }
2014: /* Overwrites A. Can only handle full-rank problems with m>=n */
2015: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2016: {
2017: PetscBool debug = PETSC_FALSE;
2018: PetscBLASInt M, N, K, lda, ldb, ldwork;
2019: PetscScalar *R, *Q, *Aback, Alpha;
2021: PetscFunctionBegin;
2022: if (debug) {
2023: PetscCall(PetscMalloc1(m * n, &Aback));
2024: PetscCall(PetscArraycpy(Aback, A, m * n));
2025: }
2027: PetscCall(PetscBLASIntCast(m, &M));
2028: PetscCall(PetscBLASIntCast(n, &N));
2029: PetscCall(PetscBLASIntCast(mstride, &lda));
2030: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2031: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2032: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2033: PetscCall(PetscFPTrapPop());
2034: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2036: /* Extract an explicit representation of Q */
2037: Q = Ainv;
2038: PetscCall(PetscArraycpy(Q, A, mstride * n));
2039: K = N; /* full rank */
2040: PetscCallLAPACKInfo("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2042: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2043: Alpha = 1.0;
2044: ldb = lda;
2045: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2046: /* Ainv is Q, overwritten with inverse */
2048: if (debug) { /* Check that pseudo-inverse worked */
2049: PetscScalar Beta = 0.0;
2050: PetscBLASInt ldc;
2051: K = N;
2052: ldc = N;
2053: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2054: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2055: PetscCall(PetscFree(Aback));
2056: }
2057: PetscFunctionReturn(PETSC_SUCCESS);
2058: }
2060: /* Overwrites A. Can handle degenerate problems and m<n. */
2061: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2062: {
2063: PetscScalar *Brhs;
2064: PetscScalar *tmpwork;
2065: PetscReal rcond;
2066: #if PetscDefined(USE_COMPLEX)
2067: PetscInt rworkSize;
2068: PetscReal *rwork, *rtau;
2069: #endif
2070: PetscInt i, j, maxmn;
2071: PetscBLASInt M, N, lda, ldb, ldwork;
2072: PetscBLASInt nrhs, irank;
2074: PetscFunctionBegin;
2075: /* initialize to identity */
2076: tmpwork = work;
2077: Brhs = Ainv;
2078: maxmn = PetscMax(m, n);
2079: for (j = 0; j < maxmn; j++) {
2080: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2081: }
2083: PetscCall(PetscBLASIntCast(m, &M));
2084: PetscCall(PetscBLASIntCast(n, &N));
2085: PetscCall(PetscBLASIntCast(mstride, &lda));
2086: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2087: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2088: rcond = -1;
2089: nrhs = M;
2090: #if PetscDefined(USE_COMPLEX)
2091: rworkSize = 5 * PetscMin(M, N);
2092: PetscCall(PetscMalloc1(rworkSize, &rwork));
2093: PetscCall(PetscMalloc1(PetscMin(M, N), &rtau));
2094: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2095: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, rtau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2096: PetscCall(PetscFPTrapPop());
2097: PetscCall(PetscFree(rwork));
2098: for (i = 0; i < PetscMin(M, N); i++) tau[i] = rtau[i];
2099: PetscCall(PetscFree(rtau));
2100: #else
2101: nrhs = M;
2102: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2103: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, tau, &rcond, &irank, tmpwork, &ldwork, &info));
2104: PetscCall(PetscFPTrapPop());
2105: #endif
2106: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2107: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two collinear points");
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: #if 0
2112: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2113: {
2114: PetscReal grad[2] = {0, 0};
2115: const PetscInt *faces;
2116: PetscInt numFaces, f;
2118: PetscFunctionBegin;
2119: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2120: PetscCall(DMPlexGetCone(dm, cell, &faces));
2121: for (f = 0; f < numFaces; ++f) {
2122: const PetscInt *fcells;
2123: const CellGeom *cg1;
2124: const FaceGeom *fg;
2126: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2127: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2128: for (i = 0; i < 2; ++i) {
2129: PetscScalar du;
2131: if (fcells[i] == c) continue;
2132: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2133: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2134: grad[0] += fg->grad[!i][0] * du;
2135: grad[1] += fg->grad[!i][1] * du;
2136: }
2137: }
2138: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2139: PetscFunctionReturn(PETSC_SUCCESS);
2140: }
2141: #endif
2143: /*
2144: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2146: Input Parameters:
2147: + fvm - The `PetscFV` object
2148: . numFaces - The number of cell faces which are not constrained
2149: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2151: Level: developer
2153: .seealso: `PetscFV`, `PetscFVCreate()`
2154: */
2155: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2156: {
2157: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2158: const PetscBool useSVD = PETSC_TRUE;
2159: const PetscInt maxFaces = ls->maxFaces;
2160: PetscInt dim, f, d;
2162: PetscFunctionBegin;
2163: if (numFaces > maxFaces) {
2164: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2165: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2166: }
2167: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2168: for (f = 0; f < numFaces; ++f) {
2169: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2170: }
2171: /* Overwrites B with garbage, returns Binv in row-major format */
2172: if (useSVD) {
2173: PetscInt maxmn = PetscMax(numFaces, dim);
2174: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2175: /* Binv shaped in column-major, coldim=maxmn.*/
2176: for (f = 0; f < numFaces; ++f) {
2177: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2178: }
2179: } else {
2180: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2181: /* Binv shaped in row-major, rowdim=maxFaces.*/
2182: for (f = 0; f < numFaces; ++f) {
2183: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2184: }
2185: }
2186: PetscFunctionReturn(PETSC_SUCCESS);
2187: }
2189: /*
2190: neighborVol[f*2+0] contains the left geom
2191: neighborVol[f*2+1] contains the right geom
2192: */
2193: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2194: {
2195: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2196: void *rctx;
2197: PetscScalar *flux = fvm->fluxWork;
2198: const PetscScalar *constants;
2199: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2201: PetscFunctionBegin;
2202: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2203: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2204: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2205: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2206: PetscCall(PetscDSGetContext(prob, field, &rctx));
2207: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2208: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2209: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2210: for (f = 0; f < Nf; ++f) {
2211: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2212: for (d = 0; d < pdim; ++d) {
2213: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2214: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2215: }
2216: }
2217: PetscFunctionReturn(PETSC_SUCCESS);
2218: }
2220: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2221: {
2222: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2223: PetscInt dim, m, n, nrhs, minmn, maxmn;
2225: PetscFunctionBegin;
2227: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2228: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2229: ls->maxFaces = maxFaces;
2230: m = ls->maxFaces;
2231: n = dim;
2232: nrhs = ls->maxFaces;
2233: minmn = PetscMin(m, n);
2234: maxmn = PetscMax(m, n);
2235: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2236: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2237: PetscFunctionReturn(PETSC_SUCCESS);
2238: }
2240: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2241: {
2242: PetscFunctionBegin;
2243: fvm->ops->setfromoptions = NULL;
2244: fvm->ops->view = PetscFVView_LeastSquares;
2245: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2246: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2247: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2248: PetscFunctionReturn(PETSC_SUCCESS);
2249: }
2251: /*MC
2252: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2254: Level: intermediate
2256: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2257: M*/
2259: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2260: {
2261: PetscFV_LeastSquares *ls;
2263: PetscFunctionBegin;
2265: PetscCall(PetscNew(&ls));
2266: fvm->data = ls;
2268: ls->maxFaces = -1;
2269: ls->workSize = -1;
2270: ls->B = NULL;
2271: ls->Binv = NULL;
2272: ls->tau = NULL;
2273: ls->work = NULL;
2275: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2276: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2277: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2278: PetscFunctionReturn(PETSC_SUCCESS);
2279: }
2281: /*@
2282: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2284: Not Collective
2286: Input Parameters:
2287: + fvm - The `PetscFV` object
2288: - maxFaces - The maximum number of cell faces
2290: Level: intermediate
2292: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2293: @*/
2294: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2295: {
2296: PetscFunctionBegin;
2298: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2299: PetscFunctionReturn(PETSC_SUCCESS);
2300: }