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: /*@C
20: PetscLimiterRegister - Adds a new `PetscLimiter` implementation
22: Not Collective
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()`, `PetscLimiterRegisterDestroy()`
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: /*@C
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: if (lim->ops->destroy) {
88: PetscUseTypeMethod(lim, destroy);
89: lim->ops->destroy = NULL;
90: }
91: PetscCall((*r)(lim));
92: PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@C
97: PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.
99: Not Collective
101: Input Parameter:
102: . lim - The `PetscLimiter`
104: Output Parameter:
105: . name - The `PetscLimiterType`
107: Level: intermediate
109: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
110: @*/
111: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
112: {
113: PetscFunctionBegin;
115: PetscAssertPointer(name, 2);
116: PetscCall(PetscLimiterRegisterAll());
117: *name = ((PetscObject)lim)->type_name;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@C
122: PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database
124: Collective
126: Input Parameters:
127: + A - the `PetscLimiter` object to view
128: . obj - Optional object that provides the options prefix to use
129: - name - command line option name
131: Level: intermediate
133: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
134: @*/
135: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
136: {
137: PetscFunctionBegin;
139: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@C
144: PetscLimiterView - Views a `PetscLimiter`
146: Collective
148: Input Parameters:
149: + lim - the `PetscLimiter` object to view
150: - v - the viewer
152: Level: beginner
154: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
155: @*/
156: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
157: {
158: PetscFunctionBegin;
160: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
161: PetscTryTypeMethod(lim, view, v);
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*@
166: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
168: Collective
170: Input Parameter:
171: . lim - the `PetscLimiter` object to set options for
173: Level: intermediate
175: .seealso: `PetscLimiter`, ``PetscLimiterView()`
176: @*/
177: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
178: {
179: const char *defaultType;
180: char name[256];
181: PetscBool flg;
183: PetscFunctionBegin;
185: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
186: else defaultType = ((PetscObject)lim)->type_name;
187: PetscCall(PetscLimiterRegisterAll());
189: PetscObjectOptionsBegin((PetscObject)lim);
190: PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg));
191: if (flg) {
192: PetscCall(PetscLimiterSetType(lim, name));
193: } else if (!((PetscObject)lim)->type_name) {
194: PetscCall(PetscLimiterSetType(lim, defaultType));
195: }
196: PetscTryTypeMethod(lim, setfromoptions);
197: /* process any options handlers added with PetscObjectAddOptionsHandler() */
198: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
199: PetscOptionsEnd();
200: PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@C
205: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
207: Collective
209: Input Parameter:
210: . lim - the `PetscLimiter` object to setup
212: Level: intermediate
214: .seealso: `PetscLimiter`, ``PetscLimiterView()`, `PetscLimiterDestroy()`
215: @*/
216: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
217: {
218: PetscFunctionBegin;
220: PetscTryTypeMethod(lim, setup);
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*@
225: PetscLimiterDestroy - Destroys a `PetscLimiter` object
227: Collective
229: Input Parameter:
230: . lim - the `PetscLimiter` object to destroy
232: Level: beginner
234: .seealso: `PetscLimiter`, `PetscLimiterView()`
235: @*/
236: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
237: {
238: PetscFunctionBegin;
239: if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);
242: if (--((PetscObject)(*lim))->refct > 0) {
243: *lim = NULL;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
246: ((PetscObject)(*lim))->refct = 0;
248: PetscTryTypeMethod((*lim), destroy);
249: PetscCall(PetscHeaderDestroy(lim));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@
254: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
256: Collective
258: Input Parameter:
259: . comm - The communicator for the `PetscLimiter` object
261: Output Parameter:
262: . lim - The `PetscLimiter` object
264: Level: beginner
266: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
267: @*/
268: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
269: {
270: PetscLimiter l;
272: PetscFunctionBegin;
273: PetscAssertPointer(lim, 2);
274: PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
275: *lim = NULL;
276: PetscCall(PetscFVInitializePackage());
278: PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));
280: *lim = l;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: PetscLimiterLimit - Limit the flux
287: Input Parameters:
288: + lim - The `PetscLimiter`
289: - flim - The input field
291: Output Parameter:
292: . phi - The limited field
294: Level: beginner
296: Note:
297: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
298: .vb
299: The classical flux-limited formulation is psi(r) where
301: r = (u[0] - u[-1]) / (u[1] - u[0])
303: The second order TVD region is bounded by
305: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
307: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
308: phi(r)(r+1)/2 in which the reconstructed interface values are
310: u(v) = u[0] + phi(r) (grad u)[0] v
312: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
314: 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))
316: For a nicer symmetric formulation, rewrite in terms of
318: f = (u[0] - u[-1]) / (u[1] - u[-1])
320: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
322: phi(r) = phi(1/r)
324: becomes
326: w(f) = w(1-f).
328: The limiters below implement this final form w(f). The reference methods are
330: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
331: .ve
333: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
334: @*/
335: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
336: {
337: PetscFunctionBegin;
339: PetscAssertPointer(phi, 3);
340: PetscUseTypeMethod(lim, limit, flim, phi);
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
345: {
346: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
348: PetscFunctionBegin;
349: PetscCall(PetscFree(l));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
354: {
355: PetscViewerFormat format;
357: PetscFunctionBegin;
358: PetscCall(PetscViewerGetFormat(viewer, &format));
359: PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
364: {
365: PetscBool iascii;
367: PetscFunctionBegin;
370: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
371: if (iascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
376: {
377: PetscFunctionBegin;
378: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
383: {
384: PetscFunctionBegin;
385: lim->ops->view = PetscLimiterView_Sin;
386: lim->ops->destroy = PetscLimiterDestroy_Sin;
387: lim->ops->limit = PetscLimiterLimit_Sin;
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*MC
392: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
394: Level: intermediate
396: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
397: M*/
399: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
400: {
401: PetscLimiter_Sin *l;
403: PetscFunctionBegin;
405: PetscCall(PetscNew(&l));
406: lim->data = l;
408: PetscCall(PetscLimiterInitialize_Sin(lim));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
413: {
414: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
416: PetscFunctionBegin;
417: PetscCall(PetscFree(l));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
422: {
423: PetscViewerFormat format;
425: PetscFunctionBegin;
426: PetscCall(PetscViewerGetFormat(viewer, &format));
427: PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
432: {
433: PetscBool iascii;
435: PetscFunctionBegin;
438: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
439: if (iascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
444: {
445: PetscFunctionBegin;
446: *phi = 0.0;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
451: {
452: PetscFunctionBegin;
453: lim->ops->view = PetscLimiterView_Zero;
454: lim->ops->destroy = PetscLimiterDestroy_Zero;
455: lim->ops->limit = PetscLimiterLimit_Zero;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /*MC
460: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
462: Level: intermediate
464: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
465: M*/
467: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
468: {
469: PetscLimiter_Zero *l;
471: PetscFunctionBegin;
473: PetscCall(PetscNew(&l));
474: lim->data = l;
476: PetscCall(PetscLimiterInitialize_Zero(lim));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
481: {
482: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
484: PetscFunctionBegin;
485: PetscCall(PetscFree(l));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
490: {
491: PetscViewerFormat format;
493: PetscFunctionBegin;
494: PetscCall(PetscViewerGetFormat(viewer, &format));
495: PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
500: {
501: PetscBool iascii;
503: PetscFunctionBegin;
506: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
507: if (iascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
512: {
513: PetscFunctionBegin;
514: *phi = 1.0;
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
519: {
520: PetscFunctionBegin;
521: lim->ops->view = PetscLimiterView_None;
522: lim->ops->destroy = PetscLimiterDestroy_None;
523: lim->ops->limit = PetscLimiterLimit_None;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*MC
528: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
530: Level: intermediate
532: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
533: M*/
535: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
536: {
537: PetscLimiter_None *l;
539: PetscFunctionBegin;
541: PetscCall(PetscNew(&l));
542: lim->data = l;
544: PetscCall(PetscLimiterInitialize_None(lim));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
549: {
550: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
552: PetscFunctionBegin;
553: PetscCall(PetscFree(l));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
558: {
559: PetscViewerFormat format;
561: PetscFunctionBegin;
562: PetscCall(PetscViewerGetFormat(viewer, &format));
563: PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
568: {
569: PetscBool iascii;
571: PetscFunctionBegin;
574: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
575: if (iascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
580: {
581: PetscFunctionBegin;
582: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
587: {
588: PetscFunctionBegin;
589: lim->ops->view = PetscLimiterView_Minmod;
590: lim->ops->destroy = PetscLimiterDestroy_Minmod;
591: lim->ops->limit = PetscLimiterLimit_Minmod;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: /*MC
596: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
598: Level: intermediate
600: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
601: M*/
603: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
604: {
605: PetscLimiter_Minmod *l;
607: PetscFunctionBegin;
609: PetscCall(PetscNew(&l));
610: lim->data = l;
612: PetscCall(PetscLimiterInitialize_Minmod(lim));
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
617: {
618: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
620: PetscFunctionBegin;
621: PetscCall(PetscFree(l));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
626: {
627: PetscViewerFormat format;
629: PetscFunctionBegin;
630: PetscCall(PetscViewerGetFormat(viewer, &format));
631: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
636: {
637: PetscBool iascii;
639: PetscFunctionBegin;
642: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
643: if (iascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
648: {
649: PetscFunctionBegin;
650: *phi = PetscMax(0, 4 * f * (1 - f));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
655: {
656: PetscFunctionBegin;
657: lim->ops->view = PetscLimiterView_VanLeer;
658: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
659: lim->ops->limit = PetscLimiterLimit_VanLeer;
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*MC
664: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
666: Level: intermediate
668: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
669: M*/
671: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
672: {
673: PetscLimiter_VanLeer *l;
675: PetscFunctionBegin;
677: PetscCall(PetscNew(&l));
678: lim->data = l;
680: PetscCall(PetscLimiterInitialize_VanLeer(lim));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
685: {
686: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
688: PetscFunctionBegin;
689: PetscCall(PetscFree(l));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
694: {
695: PetscViewerFormat format;
697: PetscFunctionBegin;
698: PetscCall(PetscViewerGetFormat(viewer, &format));
699: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
704: {
705: PetscBool iascii;
707: PetscFunctionBegin;
710: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
711: if (iascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
716: {
717: PetscFunctionBegin;
718: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
723: {
724: PetscFunctionBegin;
725: lim->ops->view = PetscLimiterView_VanAlbada;
726: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
727: lim->ops->limit = PetscLimiterLimit_VanAlbada;
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*MC
732: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
734: Level: intermediate
736: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
737: M*/
739: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
740: {
741: PetscLimiter_VanAlbada *l;
743: PetscFunctionBegin;
745: PetscCall(PetscNew(&l));
746: lim->data = l;
748: PetscCall(PetscLimiterInitialize_VanAlbada(lim));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
753: {
754: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
756: PetscFunctionBegin;
757: PetscCall(PetscFree(l));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
762: {
763: PetscViewerFormat format;
765: PetscFunctionBegin;
766: PetscCall(PetscViewerGetFormat(viewer, &format));
767: PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
772: {
773: PetscBool iascii;
775: PetscFunctionBegin;
778: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
779: if (iascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
784: {
785: PetscFunctionBegin;
786: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
791: {
792: PetscFunctionBegin;
793: lim->ops->view = PetscLimiterView_Superbee;
794: lim->ops->destroy = PetscLimiterDestroy_Superbee;
795: lim->ops->limit = PetscLimiterLimit_Superbee;
796: PetscFunctionReturn(PETSC_SUCCESS);
797: }
799: /*MC
800: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
802: Level: intermediate
804: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
805: M*/
807: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
808: {
809: PetscLimiter_Superbee *l;
811: PetscFunctionBegin;
813: PetscCall(PetscNew(&l));
814: lim->data = l;
816: PetscCall(PetscLimiterInitialize_Superbee(lim));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
821: {
822: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
824: PetscFunctionBegin;
825: PetscCall(PetscFree(l));
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
830: {
831: PetscViewerFormat format;
833: PetscFunctionBegin;
834: PetscCall(PetscViewerGetFormat(viewer, &format));
835: PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
840: {
841: PetscBool iascii;
843: PetscFunctionBegin;
846: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
847: if (iascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /* aka Barth-Jespersen */
852: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
853: {
854: PetscFunctionBegin;
855: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
860: {
861: PetscFunctionBegin;
862: lim->ops->view = PetscLimiterView_MC;
863: lim->ops->destroy = PetscLimiterDestroy_MC;
864: lim->ops->limit = PetscLimiterLimit_MC;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*MC
869: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
871: Level: intermediate
873: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
874: M*/
876: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
877: {
878: PetscLimiter_MC *l;
880: PetscFunctionBegin;
882: PetscCall(PetscNew(&l));
883: lim->data = l;
885: PetscCall(PetscLimiterInitialize_MC(lim));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: PetscClassId PETSCFV_CLASSID = 0;
891: PetscFunctionList PetscFVList = NULL;
892: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
894: /*@C
895: PetscFVRegister - Adds a new `PetscFV` implementation
897: Not Collective
899: Input Parameters:
900: + sname - The name of a new user-defined creation routine
901: - function - The creation routine itself
903: Example Usage:
904: .vb
905: PetscFVRegister("my_fv", MyPetscFVCreate);
906: .ve
908: Then, your PetscFV type can be chosen with the procedural interface via
909: .vb
910: PetscFVCreate(MPI_Comm, PetscFV *);
911: PetscFVSetType(PetscFV, "my_fv");
912: .ve
913: or at runtime via the option
914: .vb
915: -petscfv_type my_fv
916: .ve
918: Level: advanced
920: Note:
921: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
923: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
924: @*/
925: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
926: {
927: PetscFunctionBegin;
928: PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
929: PetscFunctionReturn(PETSC_SUCCESS);
930: }
932: /*@C
933: PetscFVSetType - Builds a particular `PetscFV`
935: Collective
937: Input Parameters:
938: + fvm - The `PetscFV` object
939: - name - The type of FVM space
941: Options Database Key:
942: . -petscfv_type <type> - Sets the `PetscFVType`; use -help for a list of available types
944: Level: intermediate
946: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
947: @*/
948: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
949: {
950: PetscErrorCode (*r)(PetscFV);
951: PetscBool match;
953: PetscFunctionBegin;
955: PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
956: if (match) PetscFunctionReturn(PETSC_SUCCESS);
958: PetscCall(PetscFVRegisterAll());
959: PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
960: PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
962: PetscTryTypeMethod(fvm, destroy);
963: fvm->ops->destroy = NULL;
965: PetscCall((*r)(fvm));
966: PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: /*@C
971: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
973: Not Collective
975: Input Parameter:
976: . fvm - The `PetscFV`
978: Output Parameter:
979: . name - The `PetscFVType` name
981: Level: intermediate
983: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
984: @*/
985: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
986: {
987: PetscFunctionBegin;
989: PetscAssertPointer(name, 2);
990: PetscCall(PetscFVRegisterAll());
991: *name = ((PetscObject)fvm)->type_name;
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: /*@C
996: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
998: Collective
1000: Input Parameters:
1001: + A - the `PetscFV` object
1002: . obj - Optional object that provides the options prefix
1003: - name - command line option name
1005: Level: intermediate
1007: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
1008: @*/
1009: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1010: {
1011: PetscFunctionBegin;
1013: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: /*@C
1018: PetscFVView - Views a `PetscFV`
1020: Collective
1022: Input Parameters:
1023: + fvm - the `PetscFV` object to view
1024: - v - the viewer
1026: Level: beginner
1028: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1029: @*/
1030: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1031: {
1032: PetscFunctionBegin;
1034: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1035: PetscTryTypeMethod(fvm, view, v);
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: /*@
1040: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
1042: Collective
1044: Input Parameter:
1045: . fvm - the `PetscFV` object to set options for
1047: Options Database Key:
1048: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1050: Level: intermediate
1052: .seealso: `PetscFV`, `PetscFVView()`
1053: @*/
1054: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1055: {
1056: const char *defaultType;
1057: char name[256];
1058: PetscBool flg;
1060: PetscFunctionBegin;
1062: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1063: else defaultType = ((PetscObject)fvm)->type_name;
1064: PetscCall(PetscFVRegisterAll());
1066: PetscObjectOptionsBegin((PetscObject)fvm);
1067: PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg));
1068: if (flg) {
1069: PetscCall(PetscFVSetType(fvm, name));
1070: } else if (!((PetscObject)fvm)->type_name) {
1071: PetscCall(PetscFVSetType(fvm, defaultType));
1072: }
1073: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1074: PetscTryTypeMethod(fvm, setfromoptions);
1075: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1076: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1077: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1078: PetscOptionsEnd();
1079: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1086: Collective
1088: Input Parameter:
1089: . fvm - the `PetscFV` object to setup
1091: Level: intermediate
1093: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1094: @*/
1095: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1096: {
1097: PetscFunctionBegin;
1099: PetscCall(PetscLimiterSetUp(fvm->limiter));
1100: PetscTryTypeMethod(fvm, setup);
1101: PetscFunctionReturn(PETSC_SUCCESS);
1102: }
1104: /*@
1105: PetscFVDestroy - Destroys a `PetscFV` object
1107: Collective
1109: Input Parameter:
1110: . fvm - the `PetscFV` object to destroy
1112: Level: beginner
1114: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1115: @*/
1116: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1117: {
1118: PetscInt i;
1120: PetscFunctionBegin;
1121: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1124: if (--((PetscObject)(*fvm))->refct > 0) {
1125: *fvm = NULL;
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1128: ((PetscObject)(*fvm))->refct = 0;
1130: for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1131: PetscCall(PetscFree((*fvm)->componentNames));
1132: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1133: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1134: PetscCall(PetscFree((*fvm)->fluxWork));
1135: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1136: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1138: PetscTryTypeMethod((*fvm), destroy);
1139: PetscCall(PetscHeaderDestroy(fvm));
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: /*@
1144: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1146: Collective
1148: Input Parameter:
1149: . comm - The communicator for the `PetscFV` object
1151: Output Parameter:
1152: . fvm - The `PetscFV` object
1154: Level: beginner
1156: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1157: @*/
1158: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1159: {
1160: PetscFV f;
1162: PetscFunctionBegin;
1163: PetscAssertPointer(fvm, 2);
1164: *fvm = NULL;
1165: PetscCall(PetscFVInitializePackage());
1167: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1168: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1170: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1171: f->numComponents = 1;
1172: f->dim = 0;
1173: f->computeGradients = PETSC_FALSE;
1174: f->fluxWork = NULL;
1175: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1177: *fvm = f;
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: /*@
1182: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1184: Logically Collective
1186: Input Parameters:
1187: + fvm - the `PetscFV` object
1188: - lim - The `PetscLimiter`
1190: Level: intermediate
1192: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1193: @*/
1194: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1195: {
1196: PetscFunctionBegin;
1199: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1200: PetscCall(PetscObjectReference((PetscObject)lim));
1201: fvm->limiter = lim;
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*@
1206: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1208: Not Collective
1210: Input Parameter:
1211: . fvm - the `PetscFV` object
1213: Output Parameter:
1214: . lim - The `PetscLimiter`
1216: Level: intermediate
1218: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1219: @*/
1220: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1221: {
1222: PetscFunctionBegin;
1224: PetscAssertPointer(lim, 2);
1225: *lim = fvm->limiter;
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: /*@
1230: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1232: Logically Collective
1234: Input Parameters:
1235: + fvm - the `PetscFV` object
1236: - comp - The number of components
1238: Level: intermediate
1240: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1241: @*/
1242: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1243: {
1244: PetscFunctionBegin;
1246: if (fvm->numComponents != comp) {
1247: PetscInt i;
1249: for (i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1250: PetscCall(PetscFree(fvm->componentNames));
1251: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1252: }
1253: fvm->numComponents = comp;
1254: PetscCall(PetscFree(fvm->fluxWork));
1255: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1256: PetscFunctionReturn(PETSC_SUCCESS);
1257: }
1259: /*@
1260: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1262: Not Collective
1264: Input Parameter:
1265: . fvm - the `PetscFV` object
1267: Output Parameter:
1268: . comp - The number of components
1270: Level: intermediate
1272: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1273: @*/
1274: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1275: {
1276: PetscFunctionBegin;
1278: PetscAssertPointer(comp, 2);
1279: *comp = fvm->numComponents;
1280: PetscFunctionReturn(PETSC_SUCCESS);
1281: }
1283: /*@C
1284: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1286: Logically Collective
1288: Input Parameters:
1289: + fvm - the `PetscFV` object
1290: . comp - the component number
1291: - name - the component name
1293: Level: intermediate
1295: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1296: @*/
1297: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1298: {
1299: PetscFunctionBegin;
1300: PetscCall(PetscFree(fvm->componentNames[comp]));
1301: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1302: PetscFunctionReturn(PETSC_SUCCESS);
1303: }
1305: /*@C
1306: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1308: Logically Collective
1310: Input Parameters:
1311: + fvm - the `PetscFV` object
1312: - comp - the component number
1314: Output Parameter:
1315: . name - the component name
1317: Level: intermediate
1319: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1320: @*/
1321: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1322: {
1323: PetscFunctionBegin;
1324: *name = fvm->componentNames[comp];
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: /*@
1329: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1331: Logically Collective
1333: Input Parameters:
1334: + fvm - the `PetscFV` object
1335: - dim - The spatial dimension
1337: Level: intermediate
1339: .seealso: `PetscFV`, ``PetscFVGetSpatialDimension()`
1340: @*/
1341: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1342: {
1343: PetscFunctionBegin;
1345: fvm->dim = dim;
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: /*@
1350: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1352: Not Collective
1354: Input Parameter:
1355: . fvm - the `PetscFV` object
1357: Output Parameter:
1358: . dim - The spatial dimension
1360: Level: intermediate
1362: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1363: @*/
1364: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1365: {
1366: PetscFunctionBegin;
1368: PetscAssertPointer(dim, 2);
1369: *dim = fvm->dim;
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: /*@
1374: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1376: Logically Collective
1378: Input Parameters:
1379: + fvm - the `PetscFV` object
1380: - computeGradients - Flag to compute cell gradients
1382: Level: intermediate
1384: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1385: @*/
1386: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1387: {
1388: PetscFunctionBegin;
1390: fvm->computeGradients = computeGradients;
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: /*@
1395: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1397: Not Collective
1399: Input Parameter:
1400: . fvm - the `PetscFV` object
1402: Output Parameter:
1403: . computeGradients - Flag to compute cell gradients
1405: Level: intermediate
1407: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1408: @*/
1409: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1410: {
1411: PetscFunctionBegin;
1413: PetscAssertPointer(computeGradients, 2);
1414: *computeGradients = fvm->computeGradients;
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: /*@
1419: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1421: Logically Collective
1423: Input Parameters:
1424: + fvm - the `PetscFV` object
1425: - q - The `PetscQuadrature`
1427: Level: intermediate
1429: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1430: @*/
1431: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1432: {
1433: PetscFunctionBegin;
1435: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1436: PetscCall(PetscObjectReference((PetscObject)q));
1437: fvm->quadrature = q;
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: /*@
1442: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1444: Not Collective
1446: Input Parameter:
1447: . fvm - the `PetscFV` object
1449: Output Parameter:
1450: . q - The `PetscQuadrature`
1452: Level: intermediate
1454: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1455: @*/
1456: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1457: {
1458: PetscFunctionBegin;
1460: PetscAssertPointer(q, 2);
1461: if (!fvm->quadrature) {
1462: /* Create default 1-point quadrature */
1463: PetscReal *points, *weights;
1465: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1466: PetscCall(PetscCalloc1(fvm->dim, &points));
1467: PetscCall(PetscMalloc1(1, &weights));
1468: weights[0] = 1.0;
1469: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1470: }
1471: *q = fvm->quadrature;
1472: PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1475: /*@
1476: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1478: Not Collective
1480: Input Parameter:
1481: . fvm - The `PetscFV` object
1483: Output Parameter:
1484: . sp - The `PetscDualSpace` object
1486: Level: intermediate
1488: Developer Notes:
1489: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1491: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1492: @*/
1493: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1494: {
1495: PetscFunctionBegin;
1497: PetscAssertPointer(sp, 2);
1498: if (!fvm->dualSpace) {
1499: DM K;
1500: PetscInt dim, Nc, c;
1502: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1503: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1504: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1505: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1506: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K));
1507: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1508: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1509: PetscCall(DMDestroy(&K));
1510: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1511: /* Should we be using PetscFVGetQuadrature() here? */
1512: for (c = 0; c < Nc; ++c) {
1513: PetscQuadrature qc;
1514: PetscReal *points, *weights;
1516: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1517: PetscCall(PetscCalloc1(dim, &points));
1518: PetscCall(PetscCalloc1(Nc, &weights));
1519: weights[c] = 1.0;
1520: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1521: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1522: PetscCall(PetscQuadratureDestroy(&qc));
1523: }
1524: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1525: }
1526: *sp = fvm->dualSpace;
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: /*@
1531: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1533: Not Collective
1535: Input Parameters:
1536: + fvm - The `PetscFV` object
1537: - sp - The `PetscDualSpace` object
1539: Level: intermediate
1541: Note:
1542: A simple dual space is provided automatically, and the user typically will not need to override it.
1544: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1545: @*/
1546: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1547: {
1548: PetscFunctionBegin;
1551: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1552: fvm->dualSpace = sp;
1553: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1554: PetscFunctionReturn(PETSC_SUCCESS);
1555: }
1557: /*@C
1558: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1560: Not Collective
1562: Input Parameter:
1563: . fvm - The `PetscFV` object
1565: Output Parameter:
1566: . T - The basis function values and derivatives at quadrature points
1568: Level: intermediate
1570: Note:
1571: .vb
1572: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1573: 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
1574: 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
1575: .ve
1577: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1578: @*/
1579: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1580: {
1581: PetscInt npoints;
1582: const PetscReal *points;
1584: PetscFunctionBegin;
1586: PetscAssertPointer(T, 2);
1587: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1588: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1589: *T = fvm->T;
1590: PetscFunctionReturn(PETSC_SUCCESS);
1591: }
1593: /*@C
1594: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1596: Not Collective
1598: Input Parameters:
1599: + fvm - The `PetscFV` object
1600: . nrepl - The number of replicas
1601: . npoints - The number of tabulation points in a replica
1602: . points - The tabulation point coordinates
1603: - K - The order of derivative to tabulate
1605: Output Parameter:
1606: . T - The basis function values and derivative at tabulation points
1608: Level: intermediate
1610: Note:
1611: .vb
1612: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1613: 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
1614: 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
1615: .ve
1617: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1618: @*/
1619: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1620: {
1621: PetscInt pdim = 1; /* Dimension of approximation space P */
1622: PetscInt cdim; /* Spatial dimension */
1623: PetscInt Nc; /* Field components */
1624: PetscInt k, p, d, c, e;
1626: PetscFunctionBegin;
1627: if (!npoints || K < 0) {
1628: *T = NULL;
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: PetscAssertPointer(points, 4);
1633: PetscAssertPointer(T, 6);
1634: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1635: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1636: PetscCall(PetscMalloc1(1, T));
1637: (*T)->K = !cdim ? 0 : K;
1638: (*T)->Nr = nrepl;
1639: (*T)->Np = npoints;
1640: (*T)->Nb = pdim;
1641: (*T)->Nc = Nc;
1642: (*T)->cdim = cdim;
1643: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1644: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1645: if (K >= 0) {
1646: for (p = 0; p < nrepl * npoints; ++p)
1647: for (d = 0; d < pdim; ++d)
1648: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1649: }
1650: if (K >= 1) {
1651: for (p = 0; p < nrepl * npoints; ++p)
1652: for (d = 0; d < pdim; ++d)
1653: for (c = 0; c < Nc; ++c)
1654: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1655: }
1656: if (K >= 2) {
1657: for (p = 0; p < nrepl * npoints; ++p)
1658: for (d = 0; d < pdim; ++d)
1659: for (c = 0; c < Nc; ++c)
1660: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1661: }
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*@C
1666: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1668: Input Parameters:
1669: + fvm - The `PetscFV` object
1670: . numFaces - The number of cell faces which are not constrained
1671: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1673: Output Parameter:
1674: . grad - the gradient
1676: Level: advanced
1678: .seealso: `PetscFV`, `PetscFVCreate()`
1679: @*/
1680: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1681: {
1682: PetscFunctionBegin;
1684: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@C
1689: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1691: Not Collective
1693: Input Parameters:
1694: + fvm - The `PetscFV` object for the field being integrated
1695: . prob - The `PetscDS` specifying the discretizations and continuum functions
1696: . field - The field being integrated
1697: . Nf - The number of faces in the chunk
1698: . fgeom - The face geometry for each face in the chunk
1699: . neighborVol - The volume for each pair of cells in the chunk
1700: . uL - The state from the cell on the left
1701: - uR - The state from the cell on the right
1703: Output Parameters:
1704: + fluxL - the left fluxes for each face
1705: - fluxR - the right fluxes for each face
1707: Level: developer
1709: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1710: @*/
1711: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1712: {
1713: PetscFunctionBegin;
1715: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1716: PetscFunctionReturn(PETSC_SUCCESS);
1717: }
1719: /*@
1720: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1721: smaller copies.
1723: Input Parameter:
1724: . fv - The initial `PetscFV`
1726: Output Parameter:
1727: . fvRef - The refined `PetscFV`
1729: Level: advanced
1731: Notes:
1732: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1733: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1734: interpolation between regularly refined meshes.
1736: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1737: @*/
1738: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1739: {
1740: PetscDualSpace Q, Qref;
1741: DM K, Kref;
1742: PetscQuadrature q, qref;
1743: DMPolytopeType ct;
1744: DMPlexTransform tr;
1745: PetscReal *v0;
1746: PetscReal *jac, *invjac;
1747: PetscInt numComp, numSubelements, s;
1749: PetscFunctionBegin;
1750: PetscCall(PetscFVGetDualSpace(fv, &Q));
1751: PetscCall(PetscFVGetQuadrature(fv, &q));
1752: PetscCall(PetscDualSpaceGetDM(Q, &K));
1753: /* Create dual space */
1754: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1755: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1756: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1757: PetscCall(DMDestroy(&Kref));
1758: PetscCall(PetscDualSpaceSetUp(Qref));
1759: /* Create volume */
1760: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1761: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1762: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1763: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1764: PetscCall(PetscFVSetUp(*fvRef));
1765: /* Create quadrature */
1766: PetscCall(DMPlexGetCellType(K, 0, &ct));
1767: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1768: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1769: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1770: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1771: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1772: for (s = 0; s < numSubelements; ++s) {
1773: PetscQuadrature qs;
1774: const PetscReal *points, *weights;
1775: PetscReal *p, *w;
1776: PetscInt dim, Nc, npoints, np;
1778: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1779: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1780: np = npoints / numSubelements;
1781: PetscCall(PetscMalloc1(np * dim, &p));
1782: PetscCall(PetscMalloc1(np * Nc, &w));
1783: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1784: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1785: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1786: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1787: PetscCall(PetscQuadratureDestroy(&qs));
1788: }
1789: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1790: PetscCall(DMPlexTransformDestroy(&tr));
1791: PetscCall(PetscQuadratureDestroy(&qref));
1792: PetscCall(PetscDualSpaceDestroy(&Qref));
1793: PetscFunctionReturn(PETSC_SUCCESS);
1794: }
1796: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1797: {
1798: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1800: PetscFunctionBegin;
1801: PetscCall(PetscFree(b));
1802: PetscFunctionReturn(PETSC_SUCCESS);
1803: }
1805: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1806: {
1807: PetscInt Nc, c;
1808: PetscViewerFormat format;
1810: PetscFunctionBegin;
1811: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1812: PetscCall(PetscViewerGetFormat(viewer, &format));
1813: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1814: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1815: for (c = 0; c < Nc; c++) {
1816: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1817: }
1818: PetscFunctionReturn(PETSC_SUCCESS);
1819: }
1821: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1822: {
1823: PetscBool iascii;
1825: PetscFunctionBegin;
1828: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1829: if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1830: PetscFunctionReturn(PETSC_SUCCESS);
1831: }
1833: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1834: {
1835: PetscFunctionBegin;
1836: PetscFunctionReturn(PETSC_SUCCESS);
1837: }
1839: /*
1840: neighborVol[f*2+0] contains the left geom
1841: neighborVol[f*2+1] contains the right geom
1842: */
1843: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1844: {
1845: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1846: void *rctx;
1847: PetscScalar *flux = fvm->fluxWork;
1848: const PetscScalar *constants;
1849: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1851: PetscFunctionBegin;
1852: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1853: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1854: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1855: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1856: PetscCall(PetscDSGetContext(prob, field, &rctx));
1857: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1858: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1859: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1860: for (f = 0; f < Nf; ++f) {
1861: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1862: for (d = 0; d < pdim; ++d) {
1863: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1864: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1865: }
1866: }
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1871: {
1872: PetscFunctionBegin;
1873: fvm->ops->setfromoptions = NULL;
1874: fvm->ops->setup = PetscFVSetUp_Upwind;
1875: fvm->ops->view = PetscFVView_Upwind;
1876: fvm->ops->destroy = PetscFVDestroy_Upwind;
1877: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1878: PetscFunctionReturn(PETSC_SUCCESS);
1879: }
1881: /*MC
1882: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1884: Level: intermediate
1886: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1887: M*/
1889: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1890: {
1891: PetscFV_Upwind *b;
1893: PetscFunctionBegin;
1895: PetscCall(PetscNew(&b));
1896: fvm->data = b;
1898: PetscCall(PetscFVInitialize_Upwind(fvm));
1899: PetscFunctionReturn(PETSC_SUCCESS);
1900: }
1902: #include <petscblaslapack.h>
1904: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1905: {
1906: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1908: PetscFunctionBegin;
1909: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1910: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1911: PetscCall(PetscFree(ls));
1912: PetscFunctionReturn(PETSC_SUCCESS);
1913: }
1915: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1916: {
1917: PetscInt Nc, c;
1918: PetscViewerFormat format;
1920: PetscFunctionBegin;
1921: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1922: PetscCall(PetscViewerGetFormat(viewer, &format));
1923: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1924: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1925: for (c = 0; c < Nc; c++) {
1926: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1927: }
1928: PetscFunctionReturn(PETSC_SUCCESS);
1929: }
1931: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1932: {
1933: PetscBool iascii;
1935: PetscFunctionBegin;
1938: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1939: if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1944: {
1945: PetscFunctionBegin;
1946: PetscFunctionReturn(PETSC_SUCCESS);
1947: }
1949: /* Overwrites A. Can only handle full-rank problems with m>=n */
1950: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1951: {
1952: PetscBool debug = PETSC_FALSE;
1953: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
1954: PetscScalar *R, *Q, *Aback, Alpha;
1956: PetscFunctionBegin;
1957: if (debug) {
1958: PetscCall(PetscMalloc1(m * n, &Aback));
1959: PetscCall(PetscArraycpy(Aback, A, m * n));
1960: }
1962: PetscCall(PetscBLASIntCast(m, &M));
1963: PetscCall(PetscBLASIntCast(n, &N));
1964: PetscCall(PetscBLASIntCast(mstride, &lda));
1965: PetscCall(PetscBLASIntCast(worksize, &ldwork));
1966: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1967: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1968: PetscCall(PetscFPTrapPop());
1969: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
1970: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1972: /* Extract an explicit representation of Q */
1973: Q = Ainv;
1974: PetscCall(PetscArraycpy(Q, A, mstride * n));
1975: K = N; /* full rank */
1976: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
1977: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
1979: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1980: Alpha = 1.0;
1981: ldb = lda;
1982: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1983: /* Ainv is Q, overwritten with inverse */
1985: if (debug) { /* Check that pseudo-inverse worked */
1986: PetscScalar Beta = 0.0;
1987: PetscBLASInt ldc;
1988: K = N;
1989: ldc = N;
1990: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1991: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
1992: PetscCall(PetscFree(Aback));
1993: }
1994: PetscFunctionReturn(PETSC_SUCCESS);
1995: }
1997: /* Overwrites A. Can handle degenerate problems and m<n. */
1998: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1999: {
2000: PetscScalar *Brhs;
2001: PetscScalar *tmpwork;
2002: PetscReal rcond;
2003: #if defined(PETSC_USE_COMPLEX)
2004: PetscInt rworkSize;
2005: PetscReal *rwork;
2006: #endif
2007: PetscInt i, j, maxmn;
2008: PetscBLASInt M, N, lda, ldb, ldwork;
2009: PetscBLASInt nrhs, irank, info;
2011: PetscFunctionBegin;
2012: /* initialize to identity */
2013: tmpwork = work;
2014: Brhs = Ainv;
2015: maxmn = PetscMax(m, n);
2016: for (j = 0; j < maxmn; j++) {
2017: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2018: }
2020: PetscCall(PetscBLASIntCast(m, &M));
2021: PetscCall(PetscBLASIntCast(n, &N));
2022: PetscCall(PetscBLASIntCast(mstride, &lda));
2023: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2024: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2025: rcond = -1;
2026: nrhs = M;
2027: #if defined(PETSC_USE_COMPLEX)
2028: rworkSize = 5 * PetscMin(M, N);
2029: PetscCall(PetscMalloc1(rworkSize, &rwork));
2030: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2031: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2032: PetscCall(PetscFPTrapPop());
2033: PetscCall(PetscFree(rwork));
2034: #else
2035: nrhs = M;
2036: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2037: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2038: PetscCall(PetscFPTrapPop());
2039: #endif
2040: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2041: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2042: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2043: PetscFunctionReturn(PETSC_SUCCESS);
2044: }
2046: #if 0
2047: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2048: {
2049: PetscReal grad[2] = {0, 0};
2050: const PetscInt *faces;
2051: PetscInt numFaces, f;
2053: PetscFunctionBegin;
2054: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2055: PetscCall(DMPlexGetCone(dm, cell, &faces));
2056: for (f = 0; f < numFaces; ++f) {
2057: const PetscInt *fcells;
2058: const CellGeom *cg1;
2059: const FaceGeom *fg;
2061: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2062: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2063: for (i = 0; i < 2; ++i) {
2064: PetscScalar du;
2066: if (fcells[i] == c) continue;
2067: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2068: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2069: grad[0] += fg->grad[!i][0] * du;
2070: grad[1] += fg->grad[!i][1] * du;
2071: }
2072: }
2073: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2074: PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2076: #endif
2078: /*
2079: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2081: Input Parameters:
2082: + fvm - The `PetscFV` object
2083: . numFaces - The number of cell faces which are not constrained
2084: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2086: Level: developer
2088: .seealso: `PetscFV`, `PetscFVCreate()`
2089: */
2090: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2091: {
2092: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2093: const PetscBool useSVD = PETSC_TRUE;
2094: const PetscInt maxFaces = ls->maxFaces;
2095: PetscInt dim, f, d;
2097: PetscFunctionBegin;
2098: if (numFaces > maxFaces) {
2099: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2100: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2101: }
2102: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2103: for (f = 0; f < numFaces; ++f) {
2104: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2105: }
2106: /* Overwrites B with garbage, returns Binv in row-major format */
2107: if (useSVD) {
2108: PetscInt maxmn = PetscMax(numFaces, dim);
2109: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2110: /* Binv shaped in column-major, coldim=maxmn.*/
2111: for (f = 0; f < numFaces; ++f) {
2112: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2113: }
2114: } else {
2115: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2116: /* Binv shaped in row-major, rowdim=maxFaces.*/
2117: for (f = 0; f < numFaces; ++f) {
2118: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2119: }
2120: }
2121: PetscFunctionReturn(PETSC_SUCCESS);
2122: }
2124: /*
2125: neighborVol[f*2+0] contains the left geom
2126: neighborVol[f*2+1] contains the right geom
2127: */
2128: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2129: {
2130: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2131: void *rctx;
2132: PetscScalar *flux = fvm->fluxWork;
2133: const PetscScalar *constants;
2134: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2136: PetscFunctionBegin;
2137: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2138: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2139: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2140: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2141: PetscCall(PetscDSGetContext(prob, field, &rctx));
2142: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2143: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2144: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2145: for (f = 0; f < Nf; ++f) {
2146: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2147: for (d = 0; d < pdim; ++d) {
2148: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2149: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2150: }
2151: }
2152: PetscFunctionReturn(PETSC_SUCCESS);
2153: }
2155: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2156: {
2157: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2158: PetscInt dim, m, n, nrhs, minmn, maxmn;
2160: PetscFunctionBegin;
2162: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2163: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2164: ls->maxFaces = maxFaces;
2165: m = ls->maxFaces;
2166: n = dim;
2167: nrhs = ls->maxFaces;
2168: minmn = PetscMin(m, n);
2169: maxmn = PetscMax(m, n);
2170: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2171: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2172: PetscFunctionReturn(PETSC_SUCCESS);
2173: }
2175: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2176: {
2177: PetscFunctionBegin;
2178: fvm->ops->setfromoptions = NULL;
2179: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2180: fvm->ops->view = PetscFVView_LeastSquares;
2181: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2182: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2183: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2184: PetscFunctionReturn(PETSC_SUCCESS);
2185: }
2187: /*MC
2188: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2190: Level: intermediate
2192: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2193: M*/
2195: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2196: {
2197: PetscFV_LeastSquares *ls;
2199: PetscFunctionBegin;
2201: PetscCall(PetscNew(&ls));
2202: fvm->data = ls;
2204: ls->maxFaces = -1;
2205: ls->workSize = -1;
2206: ls->B = NULL;
2207: ls->Binv = NULL;
2208: ls->tau = NULL;
2209: ls->work = NULL;
2211: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2212: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2213: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: /*@
2218: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2220: Not Collective
2222: Input Parameters:
2223: + fvm - The `PetscFV` object
2224: - maxFaces - The maximum number of cell faces
2226: Level: intermediate
2228: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2229: @*/
2230: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2231: {
2232: PetscFunctionBegin;
2234: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2235: PetscFunctionReturn(PETSC_SUCCESS);
2236: }