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: Sample 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;
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;
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;
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: Sample 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;
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: `PetscFVSet`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1157: @*/
1158: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1159: {
1160: PetscFV f;
1162: PetscFunctionBegin;
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;
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;
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
1309: Input Parameters:
1310: + fvm - the `PetscFV` object
1311: - comp - the component number
1313: Output Parameter:
1314: . name - the component name
1316: Level: intermediate
1318: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1319: @*/
1320: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1321: {
1322: PetscFunctionBegin;
1323: *name = fvm->componentNames[comp];
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: /*@
1328: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1330: Logically Collective
1332: Input Parameters:
1333: + fvm - the `PetscFV` object
1334: - dim - The spatial dimension
1336: Level: intermediate
1338: .seealso: `PetscFV`, ``PetscFVGetSpatialDimension()`
1339: @*/
1340: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1341: {
1342: PetscFunctionBegin;
1344: fvm->dim = dim;
1345: PetscFunctionReturn(PETSC_SUCCESS);
1346: }
1348: /*@
1349: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1351: Not Collective
1353: Input Parameter:
1354: . fvm - the `PetscFV` object
1356: Output Parameter:
1357: . dim - The spatial dimension
1359: Level: intermediate
1361: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1362: @*/
1363: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1364: {
1365: PetscFunctionBegin;
1368: *dim = fvm->dim;
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }
1372: /*@
1373: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1375: Logically Collective
1377: Input Parameters:
1378: + fvm - the `PetscFV` object
1379: - computeGradients - Flag to compute cell gradients
1381: Level: intermediate
1383: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1384: @*/
1385: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1386: {
1387: PetscFunctionBegin;
1389: fvm->computeGradients = computeGradients;
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: /*@
1394: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1396: Not Collective
1398: Input Parameter:
1399: . fvm - the `PetscFV` object
1401: Output Parameter:
1402: . computeGradients - Flag to compute cell gradients
1404: Level: intermediate
1406: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1407: @*/
1408: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1409: {
1410: PetscFunctionBegin;
1413: *computeGradients = fvm->computeGradients;
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@
1418: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1420: Logically Collective
1422: Input Parameters:
1423: + fvm - the `PetscFV` object
1424: - q - The `PetscQuadrature`
1426: Level: intermediate
1428: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1429: @*/
1430: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1431: {
1432: PetscFunctionBegin;
1434: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1435: PetscCall(PetscObjectReference((PetscObject)q));
1436: fvm->quadrature = q;
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: /*@
1441: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1443: Not Collective
1445: Input Parameter:
1446: . fvm - the `PetscFV` object
1448: Output Parameter:
1449: . lim - The `PetscQuadrature`
1451: Level: intermediate
1453: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1454: @*/
1455: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1456: {
1457: PetscFunctionBegin;
1460: if (!fvm->quadrature) {
1461: /* Create default 1-point quadrature */
1462: PetscReal *points, *weights;
1464: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1465: PetscCall(PetscCalloc1(fvm->dim, &points));
1466: PetscCall(PetscMalloc1(1, &weights));
1467: weights[0] = 1.0;
1468: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1469: }
1470: *q = fvm->quadrature;
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: /*@
1475: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1477: Not Collective
1479: Input Parameter:
1480: . fvm - The `PetscFV` object
1482: Output Parameter:
1483: . sp - The `PetscDualSpace` object
1485: Level: intermediate
1487: Developer Note:
1488: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1490: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1491: @*/
1492: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1493: {
1494: PetscFunctionBegin;
1497: if (!fvm->dualSpace) {
1498: DM K;
1499: PetscInt dim, Nc, c;
1501: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1502: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1503: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1504: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1505: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K));
1506: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1507: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1508: PetscCall(DMDestroy(&K));
1509: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1510: /* Should we be using PetscFVGetQuadrature() here? */
1511: for (c = 0; c < Nc; ++c) {
1512: PetscQuadrature qc;
1513: PetscReal *points, *weights;
1515: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1516: PetscCall(PetscCalloc1(dim, &points));
1517: PetscCall(PetscCalloc1(Nc, &weights));
1518: weights[c] = 1.0;
1519: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1520: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1521: PetscCall(PetscQuadratureDestroy(&qc));
1522: }
1523: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1524: }
1525: *sp = fvm->dualSpace;
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1532: Not Collective
1534: Input Parameters:
1535: + fvm - The `PetscFV` object
1536: - sp - The `PetscDualSpace` object
1538: Level: intermediate
1540: Note:
1541: A simple dual space is provided automatically, and the user typically will not need to override it.
1543: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1544: @*/
1545: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1546: {
1547: PetscFunctionBegin;
1550: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1551: fvm->dualSpace = sp;
1552: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: /*@C
1557: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1559: Not Collective
1561: Input Parameter:
1562: . fvm - The `PetscFV` object
1564: Output Parameter:
1565: . T - The basis function values and derivatives at quadrature points
1567: Level: intermediate
1569: Note:
1570: .vb
1571: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1572: 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
1573: 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
1574: .ve
1576: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1577: @*/
1578: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1579: {
1580: PetscInt npoints;
1581: const PetscReal *points;
1583: PetscFunctionBegin;
1586: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1587: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1588: *T = fvm->T;
1589: PetscFunctionReturn(PETSC_SUCCESS);
1590: }
1592: /*@C
1593: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1595: Not Collective
1597: Input Parameters:
1598: + fvm - The `PetscFV` object
1599: . nrepl - The number of replicas
1600: . npoints - The number of tabulation points in a replica
1601: . points - The tabulation point coordinates
1602: - K - The order of derivative to tabulate
1604: Output Parameter:
1605: . T - The basis function values and derivative at tabulation points
1607: Level: intermediate
1609: Note:
1610: .vb
1611: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1612: 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
1613: 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
1614: .ve
1616: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1617: @*/
1618: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1619: {
1620: PetscInt pdim = 1; /* Dimension of approximation space P */
1621: PetscInt cdim; /* Spatial dimension */
1622: PetscInt Nc; /* Field components */
1623: PetscInt k, p, d, c, e;
1625: PetscFunctionBegin;
1626: if (!npoints || K < 0) {
1627: *T = NULL;
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1633: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1634: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1635: PetscCall(PetscMalloc1(1, T));
1636: (*T)->K = !cdim ? 0 : K;
1637: (*T)->Nr = nrepl;
1638: (*T)->Np = npoints;
1639: (*T)->Nb = pdim;
1640: (*T)->Nc = Nc;
1641: (*T)->cdim = cdim;
1642: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1643: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1644: if (K >= 0) {
1645: for (p = 0; p < nrepl * npoints; ++p)
1646: for (d = 0; d < pdim; ++d)
1647: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1648: }
1649: if (K >= 1) {
1650: for (p = 0; p < nrepl * npoints; ++p)
1651: for (d = 0; d < pdim; ++d)
1652: for (c = 0; c < Nc; ++c)
1653: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1654: }
1655: if (K >= 2) {
1656: for (p = 0; p < nrepl * npoints; ++p)
1657: for (d = 0; d < pdim; ++d)
1658: for (c = 0; c < Nc; ++c)
1659: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1660: }
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: /*@C
1665: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1667: Input Parameters:
1668: + fvm - The `PetscFV` object
1669: . numFaces - The number of cell faces which are not constrained
1670: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1672: Level: advanced
1674: .seealso: `PetscFV`, `PetscFVCreate()`
1675: @*/
1676: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1677: {
1678: PetscFunctionBegin;
1680: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: /*@C
1685: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1687: Not Collective
1689: Input Parameters:
1690: + fvm - The `PetscFV` object for the field being integrated
1691: . prob - The `PetscDS` specifying the discretizations and continuum functions
1692: . field - The field being integrated
1693: . Nf - The number of faces in the chunk
1694: . fgeom - The face geometry for each face in the chunk
1695: . neighborVol - The volume for each pair of cells in the chunk
1696: . uL - The state from the cell on the left
1697: - uR - The state from the cell on the right
1699: Output Parameters:
1700: + fluxL - the left fluxes for each face
1701: - fluxR - the right fluxes for each face
1703: Level: developer
1705: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1706: @*/
1707: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1708: {
1709: PetscFunctionBegin;
1711: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1712: PetscFunctionReturn(PETSC_SUCCESS);
1713: }
1715: /*@
1716: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into smaller copies. This is typically used
1717: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1718: sparsity). It is also used to create an interpolation between regularly refined meshes.
1720: Input Parameter:
1721: . fv - The initial `PetscFV`
1723: Output Parameter:
1724: . fvRef - The refined `PetscFV`
1726: Level: advanced
1728: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1729: @*/
1730: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1731: {
1732: PetscDualSpace Q, Qref;
1733: DM K, Kref;
1734: PetscQuadrature q, qref;
1735: DMPolytopeType ct;
1736: DMPlexTransform tr;
1737: PetscReal *v0;
1738: PetscReal *jac, *invjac;
1739: PetscInt numComp, numSubelements, s;
1741: PetscFunctionBegin;
1742: PetscCall(PetscFVGetDualSpace(fv, &Q));
1743: PetscCall(PetscFVGetQuadrature(fv, &q));
1744: PetscCall(PetscDualSpaceGetDM(Q, &K));
1745: /* Create dual space */
1746: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1747: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1748: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1749: PetscCall(DMDestroy(&Kref));
1750: PetscCall(PetscDualSpaceSetUp(Qref));
1751: /* Create volume */
1752: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1753: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1754: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1755: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1756: PetscCall(PetscFVSetUp(*fvRef));
1757: /* Create quadrature */
1758: PetscCall(DMPlexGetCellType(K, 0, &ct));
1759: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1760: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1761: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1762: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1763: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1764: for (s = 0; s < numSubelements; ++s) {
1765: PetscQuadrature qs;
1766: const PetscReal *points, *weights;
1767: PetscReal *p, *w;
1768: PetscInt dim, Nc, npoints, np;
1770: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1771: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1772: np = npoints / numSubelements;
1773: PetscCall(PetscMalloc1(np * dim, &p));
1774: PetscCall(PetscMalloc1(np * Nc, &w));
1775: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1776: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1777: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1778: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1779: PetscCall(PetscQuadratureDestroy(&qs));
1780: }
1781: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1782: PetscCall(DMPlexTransformDestroy(&tr));
1783: PetscCall(PetscQuadratureDestroy(&qref));
1784: PetscCall(PetscDualSpaceDestroy(&Qref));
1785: PetscFunctionReturn(PETSC_SUCCESS);
1786: }
1788: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1789: {
1790: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1792: PetscFunctionBegin;
1793: PetscCall(PetscFree(b));
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1798: {
1799: PetscInt Nc, c;
1800: PetscViewerFormat format;
1802: PetscFunctionBegin;
1803: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1804: PetscCall(PetscViewerGetFormat(viewer, &format));
1805: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1806: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1807: for (c = 0; c < Nc; c++) {
1808: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1809: }
1810: PetscFunctionReturn(PETSC_SUCCESS);
1811: }
1813: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1814: {
1815: PetscBool iascii;
1817: PetscFunctionBegin;
1820: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1821: if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1822: PetscFunctionReturn(PETSC_SUCCESS);
1823: }
1825: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1826: {
1827: PetscFunctionBegin;
1828: PetscFunctionReturn(PETSC_SUCCESS);
1829: }
1831: /*
1832: neighborVol[f*2+0] contains the left geom
1833: neighborVol[f*2+1] contains the right geom
1834: */
1835: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1836: {
1837: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1838: void *rctx;
1839: PetscScalar *flux = fvm->fluxWork;
1840: const PetscScalar *constants;
1841: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1843: PetscFunctionBegin;
1844: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1845: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1846: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1847: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1848: PetscCall(PetscDSGetContext(prob, field, &rctx));
1849: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1850: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1851: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1852: for (f = 0; f < Nf; ++f) {
1853: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1854: for (d = 0; d < pdim; ++d) {
1855: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1856: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1857: }
1858: }
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1863: {
1864: PetscFunctionBegin;
1865: fvm->ops->setfromoptions = NULL;
1866: fvm->ops->setup = PetscFVSetUp_Upwind;
1867: fvm->ops->view = PetscFVView_Upwind;
1868: fvm->ops->destroy = PetscFVDestroy_Upwind;
1869: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1873: /*MC
1874: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1876: Level: intermediate
1878: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1879: M*/
1881: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1882: {
1883: PetscFV_Upwind *b;
1885: PetscFunctionBegin;
1887: PetscCall(PetscNew(&b));
1888: fvm->data = b;
1890: PetscCall(PetscFVInitialize_Upwind(fvm));
1891: PetscFunctionReturn(PETSC_SUCCESS);
1892: }
1894: #include <petscblaslapack.h>
1896: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1897: {
1898: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1900: PetscFunctionBegin;
1901: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1902: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1903: PetscCall(PetscFree(ls));
1904: PetscFunctionReturn(PETSC_SUCCESS);
1905: }
1907: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1908: {
1909: PetscInt Nc, c;
1910: PetscViewerFormat format;
1912: PetscFunctionBegin;
1913: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1914: PetscCall(PetscViewerGetFormat(viewer, &format));
1915: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1916: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1917: for (c = 0; c < Nc; c++) {
1918: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1919: }
1920: PetscFunctionReturn(PETSC_SUCCESS);
1921: }
1923: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1924: {
1925: PetscBool iascii;
1927: PetscFunctionBegin;
1930: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1931: if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
1932: PetscFunctionReturn(PETSC_SUCCESS);
1933: }
1935: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1936: {
1937: PetscFunctionBegin;
1938: PetscFunctionReturn(PETSC_SUCCESS);
1939: }
1941: /* Overwrites A. Can only handle full-rank problems with m>=n */
1942: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1943: {
1944: PetscBool debug = PETSC_FALSE;
1945: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
1946: PetscScalar *R, *Q, *Aback, Alpha;
1948: PetscFunctionBegin;
1949: if (debug) {
1950: PetscCall(PetscMalloc1(m * n, &Aback));
1951: PetscCall(PetscArraycpy(Aback, A, m * n));
1952: }
1954: PetscCall(PetscBLASIntCast(m, &M));
1955: PetscCall(PetscBLASIntCast(n, &N));
1956: PetscCall(PetscBLASIntCast(mstride, &lda));
1957: PetscCall(PetscBLASIntCast(worksize, &ldwork));
1958: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1959: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1960: PetscCall(PetscFPTrapPop());
1961: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
1962: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1964: /* Extract an explicit representation of Q */
1965: Q = Ainv;
1966: PetscCall(PetscArraycpy(Q, A, mstride * n));
1967: K = N; /* full rank */
1968: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
1969: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
1971: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1972: Alpha = 1.0;
1973: ldb = lda;
1974: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1975: /* Ainv is Q, overwritten with inverse */
1977: if (debug) { /* Check that pseudo-inverse worked */
1978: PetscScalar Beta = 0.0;
1979: PetscBLASInt ldc;
1980: K = N;
1981: ldc = N;
1982: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1983: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
1984: PetscCall(PetscFree(Aback));
1985: }
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: /* Overwrites A. Can handle degenerate problems and m<n. */
1990: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1991: {
1992: PetscScalar *Brhs;
1993: PetscScalar *tmpwork;
1994: PetscReal rcond;
1995: #if defined(PETSC_USE_COMPLEX)
1996: PetscInt rworkSize;
1997: PetscReal *rwork;
1998: #endif
1999: PetscInt i, j, maxmn;
2000: PetscBLASInt M, N, lda, ldb, ldwork;
2001: PetscBLASInt nrhs, irank, info;
2003: PetscFunctionBegin;
2004: /* initialize to identity */
2005: tmpwork = work;
2006: Brhs = Ainv;
2007: maxmn = PetscMax(m, n);
2008: for (j = 0; j < maxmn; j++) {
2009: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2010: }
2012: PetscCall(PetscBLASIntCast(m, &M));
2013: PetscCall(PetscBLASIntCast(n, &N));
2014: PetscCall(PetscBLASIntCast(mstride, &lda));
2015: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2016: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2017: rcond = -1;
2018: nrhs = M;
2019: #if defined(PETSC_USE_COMPLEX)
2020: rworkSize = 5 * PetscMin(M, N);
2021: PetscCall(PetscMalloc1(rworkSize, &rwork));
2022: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2023: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2024: PetscCall(PetscFPTrapPop());
2025: PetscCall(PetscFree(rwork));
2026: #else
2027: nrhs = M;
2028: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2029: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2030: PetscCall(PetscFPTrapPop());
2031: #endif
2032: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2033: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2034: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2035: PetscFunctionReturn(PETSC_SUCCESS);
2036: }
2038: #if 0
2039: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2040: {
2041: PetscReal grad[2] = {0, 0};
2042: const PetscInt *faces;
2043: PetscInt numFaces, f;
2045: PetscFunctionBegin;
2046: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2047: PetscCall(DMPlexGetCone(dm, cell, &faces));
2048: for (f = 0; f < numFaces; ++f) {
2049: const PetscInt *fcells;
2050: const CellGeom *cg1;
2051: const FaceGeom *fg;
2053: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2054: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2055: for (i = 0; i < 2; ++i) {
2056: PetscScalar du;
2058: if (fcells[i] == c) continue;
2059: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2060: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2061: grad[0] += fg->grad[!i][0] * du;
2062: grad[1] += fg->grad[!i][1] * du;
2063: }
2064: }
2065: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2068: #endif
2070: /*
2071: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2073: Input Parameters:
2074: + fvm - The `PetscFV` object
2075: . numFaces - The number of cell faces which are not constrained
2076: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2078: Level: developer
2080: .seealso: `PetscFV`, `PetscFVCreate()`
2081: */
2082: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2083: {
2084: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2085: const PetscBool useSVD = PETSC_TRUE;
2086: const PetscInt maxFaces = ls->maxFaces;
2087: PetscInt dim, f, d;
2089: PetscFunctionBegin;
2090: if (numFaces > maxFaces) {
2091: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2092: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2093: }
2094: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2095: for (f = 0; f < numFaces; ++f) {
2096: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2097: }
2098: /* Overwrites B with garbage, returns Binv in row-major format */
2099: if (useSVD) {
2100: PetscInt maxmn = PetscMax(numFaces, dim);
2101: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2102: /* Binv shaped in column-major, coldim=maxmn.*/
2103: for (f = 0; f < numFaces; ++f) {
2104: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2105: }
2106: } else {
2107: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2108: /* Binv shaped in row-major, rowdim=maxFaces.*/
2109: for (f = 0; f < numFaces; ++f) {
2110: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2111: }
2112: }
2113: PetscFunctionReturn(PETSC_SUCCESS);
2114: }
2116: /*
2117: neighborVol[f*2+0] contains the left geom
2118: neighborVol[f*2+1] contains the right geom
2119: */
2120: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2121: {
2122: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2123: void *rctx;
2124: PetscScalar *flux = fvm->fluxWork;
2125: const PetscScalar *constants;
2126: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2128: PetscFunctionBegin;
2129: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2130: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2131: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2132: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2133: PetscCall(PetscDSGetContext(prob, field, &rctx));
2134: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2135: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2136: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2137: for (f = 0; f < Nf; ++f) {
2138: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2139: for (d = 0; d < pdim; ++d) {
2140: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2141: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2142: }
2143: }
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2148: {
2149: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2150: PetscInt dim, m, n, nrhs, minmn, maxmn;
2152: PetscFunctionBegin;
2154: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2155: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2156: ls->maxFaces = maxFaces;
2157: m = ls->maxFaces;
2158: n = dim;
2159: nrhs = ls->maxFaces;
2160: minmn = PetscMin(m, n);
2161: maxmn = PetscMax(m, n);
2162: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2163: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2164: PetscFunctionReturn(PETSC_SUCCESS);
2165: }
2167: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2168: {
2169: PetscFunctionBegin;
2170: fvm->ops->setfromoptions = NULL;
2171: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2172: fvm->ops->view = PetscFVView_LeastSquares;
2173: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2174: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2175: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2176: PetscFunctionReturn(PETSC_SUCCESS);
2177: }
2179: /*MC
2180: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2182: Level: intermediate
2184: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2185: M*/
2187: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2188: {
2189: PetscFV_LeastSquares *ls;
2191: PetscFunctionBegin;
2193: PetscCall(PetscNew(&ls));
2194: fvm->data = ls;
2196: ls->maxFaces = -1;
2197: ls->workSize = -1;
2198: ls->B = NULL;
2199: ls->Binv = NULL;
2200: ls->tau = NULL;
2201: ls->work = NULL;
2203: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2204: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2205: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2206: PetscFunctionReturn(PETSC_SUCCESS);
2207: }
2209: /*@
2210: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2212: Not Collective
2214: Input parameters:
2215: + fvm - The `PetscFV` object
2216: - maxFaces - The maximum number of cell faces
2218: Level: intermediate
2220: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2221: @*/
2222: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2223: {
2224: PetscFunctionBegin;
2226: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2227: PetscFunctionReturn(PETSC_SUCCESS);
2228: }