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: + name - The name of a new user-defined creation routine
26: - create_func - The creation routine itself
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: PetscFunctionListAdd(&PetscLimiterList, sname, function);
53: return 0;
54: }
56: /*@C
57: PetscLimiterSetType - Builds a `PetscLimiter` for a given `PetscLimiterType`
59: Collective on lim
61: Input Parameters:
62: + lim - The `PetscLimiter` object
63: - name - The kind of limiter
65: Options Database Key:
66: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
68: Level: intermediate
70: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
71: @*/
72: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
73: {
74: PetscErrorCode (*r)(PetscLimiter);
75: PetscBool match;
78: PetscObjectTypeCompare((PetscObject)lim, name, &match);
79: if (match) return 0;
81: PetscLimiterRegisterAll();
82: PetscFunctionListFind(PetscLimiterList, name, &r);
85: if (lim->ops->destroy) {
86: PetscUseTypeMethod(lim, destroy);
87: lim->ops->destroy = NULL;
88: }
89: (*r)(lim);
90: PetscObjectChangeTypeName((PetscObject)lim, name);
91: return 0;
92: }
94: /*@C
95: PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.
97: Not Collective
99: Input Parameter:
100: . lim - The `PetscLimiter`
102: Output Parameter:
103: . name - The `PetscLimiterType`
105: Level: intermediate
107: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
108: @*/
109: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
110: {
113: PetscLimiterRegisterAll();
114: *name = ((PetscObject)lim)->type_name;
115: return 0;
116: }
118: /*@C
119: PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database
121: Collective on A
123: Input Parameters:
124: + A - the `PetscLimiter` object to view
125: . obj - Optional object that provides the options prefix to use
126: - name - command line option name
128: Level: intermediate
130: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
131: @*/
132: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
133: {
135: PetscObjectViewFromOptions((PetscObject)A, obj, name);
136: return 0;
137: }
139: /*@C
140: PetscLimiterView - Views a `PetscLimiter`
142: Collective on lim
144: Input Parameters:
145: + lim - the `PetscLimiter` object to view
146: - v - the viewer
148: Level: beginner
150: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
151: @*/
152: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
153: {
155: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v);
156: PetscTryTypeMethod(lim, view, v);
157: return 0;
158: }
160: /*@
161: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
163: Collective on lim
165: Input Parameter:
166: . lim - the `PetscLimiter` object to set options for
168: Level: intermediate
170: .seealso: `PetscLimiter`, ``PetscLimiterView()`
171: @*/
172: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
173: {
174: const char *defaultType;
175: char name[256];
176: PetscBool flg;
179: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
180: else defaultType = ((PetscObject)lim)->type_name;
181: PetscLimiterRegisterAll();
183: PetscObjectOptionsBegin((PetscObject)lim);
184: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
185: if (flg) {
186: PetscLimiterSetType(lim, name);
187: } else if (!((PetscObject)lim)->type_name) {
188: PetscLimiterSetType(lim, defaultType);
189: }
190: PetscTryTypeMethod(lim, setfromoptions);
191: /* process any options handlers added with PetscObjectAddOptionsHandler() */
192: PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject);
193: PetscOptionsEnd();
194: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
195: return 0;
196: }
198: /*@C
199: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
201: Collective on lim
203: Input Parameter:
204: . lim - the `PetscLimiter` object to setup
206: Level: intermediate
208: .seealso: `PetscLimiter`, ``PetscLimiterView()`, `PetscLimiterDestroy()`
209: @*/
210: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
211: {
213: PetscTryTypeMethod(lim, setup);
214: return 0;
215: }
217: /*@
218: PetscLimiterDestroy - Destroys a `PetscLimiter` object
220: Collective on lim
222: Input Parameter:
223: . lim - the `PetscLimiter` object to destroy
225: Level: beginner
227: .seealso: `PetscLimiter`, `PetscLimiterView()`
228: @*/
229: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
230: {
231: if (!*lim) return 0;
234: if (--((PetscObject)(*lim))->refct > 0) {
235: *lim = NULL;
236: return 0;
237: }
238: ((PetscObject)(*lim))->refct = 0;
240: PetscTryTypeMethod((*lim), destroy);
241: PetscHeaderDestroy(lim);
242: return 0;
243: }
245: /*@
246: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
248: Collective
250: Input Parameter:
251: . comm - The communicator for the `PetscLimiter` object
253: Output Parameter:
254: . lim - The `PetscLimiter` object
256: Level: beginner
258: .seealso: `PetscLimiter`, PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
259: @*/
260: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
261: {
262: PetscLimiter l;
265: PetscCitationsRegister(LimiterCitation, &Limitercite);
266: *lim = NULL;
267: PetscFVInitializePackage();
269: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
271: *lim = l;
272: return 0;
273: }
275: /*@
276: PetscLimiterLimit - Limit the flux
278: Input Parameters:
279: + lim - The `PetscLimiter`
280: - flim - The input field
282: Output Parameter:
283: . phi - The limited field
285: Level: beginner
287: Note:
288: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
289: .vb
290: The classical flux-limited formulation is psi(r) where
292: r = (u[0] - u[-1]) / (u[1] - u[0])
294: The second order TVD region is bounded by
296: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
298: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
299: phi(r)(r+1)/2 in which the reconstructed interface values are
301: u(v) = u[0] + phi(r) (grad u)[0] v
303: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
305: 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))
307: For a nicer symmetric formulation, rewrite in terms of
309: f = (u[0] - u[-1]) / (u[1] - u[-1])
311: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
313: phi(r) = phi(1/r)
315: becomes
317: w(f) = w(1-f).
319: The limiters below implement this final form w(f). The reference methods are
321: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
322: .ve
324: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
325: @*/
326: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
327: {
330: PetscUseTypeMethod(lim, limit, flim, phi);
331: return 0;
332: }
334: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
335: {
336: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
338: PetscFree(l);
339: return 0;
340: }
342: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
343: {
344: PetscViewerFormat format;
346: PetscViewerGetFormat(viewer, &format);
347: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
348: return 0;
349: }
351: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
352: {
353: PetscBool iascii;
357: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
358: if (iascii) PetscLimiterView_Sin_Ascii(lim, viewer);
359: return 0;
360: }
362: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
363: {
364: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
365: return 0;
366: }
368: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
369: {
370: lim->ops->view = PetscLimiterView_Sin;
371: lim->ops->destroy = PetscLimiterDestroy_Sin;
372: lim->ops->limit = PetscLimiterLimit_Sin;
373: return 0;
374: }
376: /*MC
377: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
379: Level: intermediate
381: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
382: M*/
384: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
385: {
386: PetscLimiter_Sin *l;
389: PetscNew(&l);
390: lim->data = l;
392: PetscLimiterInitialize_Sin(lim);
393: return 0;
394: }
396: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
397: {
398: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
400: PetscFree(l);
401: return 0;
402: }
404: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
405: {
406: PetscViewerFormat format;
408: PetscViewerGetFormat(viewer, &format);
409: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
410: return 0;
411: }
413: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
414: {
415: PetscBool iascii;
419: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
420: if (iascii) PetscLimiterView_Zero_Ascii(lim, viewer);
421: return 0;
422: }
424: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
425: {
426: *phi = 0.0;
427: return 0;
428: }
430: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
431: {
432: lim->ops->view = PetscLimiterView_Zero;
433: lim->ops->destroy = PetscLimiterDestroy_Zero;
434: lim->ops->limit = PetscLimiterLimit_Zero;
435: return 0;
436: }
438: /*MC
439: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
441: Level: intermediate
443: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
444: M*/
446: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
447: {
448: PetscLimiter_Zero *l;
451: PetscNew(&l);
452: lim->data = l;
454: PetscLimiterInitialize_Zero(lim);
455: return 0;
456: }
458: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
459: {
460: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
462: PetscFree(l);
463: return 0;
464: }
466: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
467: {
468: PetscViewerFormat format;
470: PetscViewerGetFormat(viewer, &format);
471: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
472: return 0;
473: }
475: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
476: {
477: PetscBool iascii;
481: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
482: if (iascii) PetscLimiterView_None_Ascii(lim, viewer);
483: return 0;
484: }
486: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
487: {
488: *phi = 1.0;
489: return 0;
490: }
492: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
493: {
494: lim->ops->view = PetscLimiterView_None;
495: lim->ops->destroy = PetscLimiterDestroy_None;
496: lim->ops->limit = PetscLimiterLimit_None;
497: return 0;
498: }
500: /*MC
501: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
503: Level: intermediate
505: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
506: M*/
508: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
509: {
510: PetscLimiter_None *l;
513: PetscNew(&l);
514: lim->data = l;
516: PetscLimiterInitialize_None(lim);
517: return 0;
518: }
520: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
521: {
522: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
524: PetscFree(l);
525: return 0;
526: }
528: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
529: {
530: PetscViewerFormat format;
532: PetscViewerGetFormat(viewer, &format);
533: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
534: return 0;
535: }
537: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
538: {
539: PetscBool iascii;
543: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
544: if (iascii) PetscLimiterView_Minmod_Ascii(lim, viewer);
545: return 0;
546: }
548: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
549: {
550: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
551: return 0;
552: }
554: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
555: {
556: lim->ops->view = PetscLimiterView_Minmod;
557: lim->ops->destroy = PetscLimiterDestroy_Minmod;
558: lim->ops->limit = PetscLimiterLimit_Minmod;
559: return 0;
560: }
562: /*MC
563: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
565: Level: intermediate
567: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
568: M*/
570: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
571: {
572: PetscLimiter_Minmod *l;
575: PetscNew(&l);
576: lim->data = l;
578: PetscLimiterInitialize_Minmod(lim);
579: return 0;
580: }
582: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
583: {
584: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
586: PetscFree(l);
587: return 0;
588: }
590: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
591: {
592: PetscViewerFormat format;
594: PetscViewerGetFormat(viewer, &format);
595: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
596: return 0;
597: }
599: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
600: {
601: PetscBool iascii;
605: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
606: if (iascii) PetscLimiterView_VanLeer_Ascii(lim, viewer);
607: return 0;
608: }
610: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
611: {
612: *phi = PetscMax(0, 4 * f * (1 - f));
613: return 0;
614: }
616: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
617: {
618: lim->ops->view = PetscLimiterView_VanLeer;
619: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
620: lim->ops->limit = PetscLimiterLimit_VanLeer;
621: return 0;
622: }
624: /*MC
625: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
627: Level: intermediate
629: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
630: M*/
632: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
633: {
634: PetscLimiter_VanLeer *l;
637: PetscNew(&l);
638: lim->data = l;
640: PetscLimiterInitialize_VanLeer(lim);
641: return 0;
642: }
644: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
645: {
646: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
648: PetscFree(l);
649: return 0;
650: }
652: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
653: {
654: PetscViewerFormat format;
656: PetscViewerGetFormat(viewer, &format);
657: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
658: return 0;
659: }
661: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
662: {
663: PetscBool iascii;
667: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
668: if (iascii) PetscLimiterView_VanAlbada_Ascii(lim, viewer);
669: return 0;
670: }
672: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
673: {
674: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
675: return 0;
676: }
678: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
679: {
680: lim->ops->view = PetscLimiterView_VanAlbada;
681: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
682: lim->ops->limit = PetscLimiterLimit_VanAlbada;
683: return 0;
684: }
686: /*MC
687: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
689: Level: intermediate
691: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
692: M*/
694: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
695: {
696: PetscLimiter_VanAlbada *l;
699: PetscNew(&l);
700: lim->data = l;
702: PetscLimiterInitialize_VanAlbada(lim);
703: return 0;
704: }
706: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
707: {
708: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
710: PetscFree(l);
711: return 0;
712: }
714: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
715: {
716: PetscViewerFormat format;
718: PetscViewerGetFormat(viewer, &format);
719: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
720: return 0;
721: }
723: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
724: {
725: PetscBool iascii;
729: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
730: if (iascii) PetscLimiterView_Superbee_Ascii(lim, viewer);
731: return 0;
732: }
734: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
735: {
736: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
737: return 0;
738: }
740: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
741: {
742: lim->ops->view = PetscLimiterView_Superbee;
743: lim->ops->destroy = PetscLimiterDestroy_Superbee;
744: lim->ops->limit = PetscLimiterLimit_Superbee;
745: return 0;
746: }
748: /*MC
749: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
751: Level: intermediate
753: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
754: M*/
756: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
757: {
758: PetscLimiter_Superbee *l;
761: PetscNew(&l);
762: lim->data = l;
764: PetscLimiterInitialize_Superbee(lim);
765: return 0;
766: }
768: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
769: {
770: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
772: PetscFree(l);
773: return 0;
774: }
776: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
777: {
778: PetscViewerFormat format;
780: PetscViewerGetFormat(viewer, &format);
781: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
782: return 0;
783: }
785: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
786: {
787: PetscBool iascii;
791: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
792: if (iascii) PetscLimiterView_MC_Ascii(lim, viewer);
793: return 0;
794: }
796: /* aka Barth-Jespersen */
797: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
798: {
799: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
800: return 0;
801: }
803: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
804: {
805: lim->ops->view = PetscLimiterView_MC;
806: lim->ops->destroy = PetscLimiterDestroy_MC;
807: lim->ops->limit = PetscLimiterLimit_MC;
808: return 0;
809: }
811: /*MC
812: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
814: Level: intermediate
816: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
817: M*/
819: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
820: {
821: PetscLimiter_MC *l;
824: PetscNew(&l);
825: lim->data = l;
827: PetscLimiterInitialize_MC(lim);
828: return 0;
829: }
831: PetscClassId PETSCFV_CLASSID = 0;
833: PetscFunctionList PetscFVList = NULL;
834: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
836: /*@C
837: PetscFVRegister - Adds a new `PetscFV` implementation
839: Not Collective
841: Input Parameters:
842: + name - The name of a new user-defined creation routine
843: - create_func - The creation routine itself
845: Sample usage:
846: .vb
847: PetscFVRegister("my_fv", MyPetscFVCreate);
848: .ve
850: Then, your PetscFV type can be chosen with the procedural interface via
851: .vb
852: PetscFVCreate(MPI_Comm, PetscFV *);
853: PetscFVSetType(PetscFV, "my_fv");
854: .ve
855: or at runtime via the option
856: .vb
857: -petscfv_type my_fv
858: .ve
860: Level: advanced
862: Note:
863: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
865: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
866: @*/
867: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
868: {
869: PetscFunctionListAdd(&PetscFVList, sname, function);
870: return 0;
871: }
873: /*@C
874: PetscFVSetType - Builds a particular `PetscFV`
876: Collective on fvm
878: Input Parameters:
879: + fvm - The `PetscFV` object
880: - name - The type of FVM space
882: Options Database Key:
883: . -petscfv_type <type> - Sets the `PetscFVType`; use -help for a list of available types
885: Level: intermediate
887: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
888: @*/
889: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
890: {
891: PetscErrorCode (*r)(PetscFV);
892: PetscBool match;
895: PetscObjectTypeCompare((PetscObject)fvm, name, &match);
896: if (match) return 0;
898: PetscFVRegisterAll();
899: PetscFunctionListFind(PetscFVList, name, &r);
902: PetscTryTypeMethod(fvm, destroy);
903: fvm->ops->destroy = NULL;
905: (*r)(fvm);
906: PetscObjectChangeTypeName((PetscObject)fvm, name);
907: return 0;
908: }
910: /*@C
911: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
913: Not Collective
915: Input Parameter:
916: . fvm - The `PetscFV`
918: Output Parameter:
919: . name - The `PetscFVType` name
921: Level: intermediate
923: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
924: @*/
925: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
926: {
929: PetscFVRegisterAll();
930: *name = ((PetscObject)fvm)->type_name;
931: return 0;
932: }
934: /*@C
935: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
937: Collective on A
939: Input Parameters:
940: + A - the `PetscFV` object
941: . obj - Optional object that provides the options prefix
942: - name - command line option name
944: Level: intermediate
946: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
947: @*/
948: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
949: {
951: PetscObjectViewFromOptions((PetscObject)A, obj, name);
952: return 0;
953: }
955: /*@C
956: PetscFVView - Views a `PetscFV`
958: Collective on fvm
960: Input Parameters:
961: + fvm - the `PetscFV` object to view
962: - v - the viewer
964: Level: beginner
966: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
967: @*/
968: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
969: {
971: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v);
972: PetscTryTypeMethod(fvm, view, v);
973: return 0;
974: }
976: /*@
977: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
979: Collective on fvm
981: Input Parameter:
982: . fvm - the `PetscFV` object to set options for
984: Options Database Key:
985: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
987: Level: intermediate
989: .seealso: `PetscFV`, `PetscFVView()`
990: @*/
991: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
992: {
993: const char *defaultType;
994: char name[256];
995: PetscBool flg;
998: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
999: else defaultType = ((PetscObject)fvm)->type_name;
1000: PetscFVRegisterAll();
1002: PetscObjectOptionsBegin((PetscObject)fvm);
1003: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1004: if (flg) {
1005: PetscFVSetType(fvm, name);
1006: } else if (!((PetscObject)fvm)->type_name) {
1007: PetscFVSetType(fvm, defaultType);
1008: }
1009: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1010: PetscTryTypeMethod(fvm, setfromoptions);
1011: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1012: PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject);
1013: PetscLimiterSetFromOptions(fvm->limiter);
1014: PetscOptionsEnd();
1015: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1016: return 0;
1017: }
1019: /*@
1020: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1022: Collective on fvm
1024: Input Parameter:
1025: . fvm - the `PetscFV` object to setup
1027: Level: intermediate
1029: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1030: @*/
1031: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1032: {
1034: PetscLimiterSetUp(fvm->limiter);
1035: PetscTryTypeMethod(fvm, setup);
1036: return 0;
1037: }
1039: /*@
1040: PetscFVDestroy - Destroys a `PetscFV` object
1042: Collective on fvm
1044: Input Parameter:
1045: . fvm - the `PetscFV` object to destroy
1047: Level: beginner
1049: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1050: @*/
1051: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1052: {
1053: PetscInt i;
1055: if (!*fvm) return 0;
1058: if (--((PetscObject)(*fvm))->refct > 0) {
1059: *fvm = NULL;
1060: return 0;
1061: }
1062: ((PetscObject)(*fvm))->refct = 0;
1064: for (i = 0; i < (*fvm)->numComponents; i++) PetscFree((*fvm)->componentNames[i]);
1065: PetscFree((*fvm)->componentNames);
1066: PetscLimiterDestroy(&(*fvm)->limiter);
1067: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1068: PetscFree((*fvm)->fluxWork);
1069: PetscQuadratureDestroy(&(*fvm)->quadrature);
1070: PetscTabulationDestroy(&(*fvm)->T);
1072: PetscTryTypeMethod((*fvm), destroy);
1073: PetscHeaderDestroy(fvm);
1074: return 0;
1075: }
1077: /*@
1078: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1080: Collective
1082: Input Parameter:
1083: . comm - The communicator for the `PetscFV` object
1085: Output Parameter:
1086: . fvm - The `PetscFV` object
1088: Level: beginner
1090: .seealso: `PetscFVSet`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1091: @*/
1092: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1093: {
1094: PetscFV f;
1097: *fvm = NULL;
1098: PetscFVInitializePackage();
1100: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1101: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1103: PetscLimiterCreate(comm, &f->limiter);
1104: f->numComponents = 1;
1105: f->dim = 0;
1106: f->computeGradients = PETSC_FALSE;
1107: f->fluxWork = NULL;
1108: PetscCalloc1(f->numComponents, &f->componentNames);
1110: *fvm = f;
1111: return 0;
1112: }
1114: /*@
1115: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1117: Logically collective on fvm
1119: Input Parameters:
1120: + fvm - the `PetscFV` object
1121: - lim - The `PetscLimiter`
1123: Level: intermediate
1125: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1126: @*/
1127: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1128: {
1131: PetscLimiterDestroy(&fvm->limiter);
1132: PetscObjectReference((PetscObject)lim);
1133: fvm->limiter = lim;
1134: return 0;
1135: }
1137: /*@
1138: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1140: Not collective
1142: Input Parameter:
1143: . fvm - the `PetscFV` object
1145: Output Parameter:
1146: . lim - The `PetscLimiter`
1148: Level: intermediate
1150: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1151: @*/
1152: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1153: {
1156: *lim = fvm->limiter;
1157: return 0;
1158: }
1160: /*@
1161: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1163: Logically collective on fvm
1165: Input Parameters:
1166: + fvm - the `PetscFV` object
1167: - comp - The number of components
1169: Level: intermediate
1171: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1172: @*/
1173: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1174: {
1176: if (fvm->numComponents != comp) {
1177: PetscInt i;
1179: for (i = 0; i < fvm->numComponents; i++) PetscFree(fvm->componentNames[i]);
1180: PetscFree(fvm->componentNames);
1181: PetscCalloc1(comp, &fvm->componentNames);
1182: }
1183: fvm->numComponents = comp;
1184: PetscFree(fvm->fluxWork);
1185: PetscMalloc1(comp, &fvm->fluxWork);
1186: return 0;
1187: }
1189: /*@
1190: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1192: Not collective
1194: Input Parameter:
1195: . fvm - the `PetscFV` object
1197: Output Parameter:
1198: , comp - The number of components
1200: Level: intermediate
1202: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1203: @*/
1204: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1205: {
1208: *comp = fvm->numComponents;
1209: return 0;
1210: }
1212: /*@C
1213: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1215: Logically collective on fvm
1217: Input Parameters:
1218: + fvm - the `PetscFV` object
1219: . comp - the component number
1220: - name - the component name
1222: Level: intermediate
1224: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1225: @*/
1226: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1227: {
1228: PetscFree(fvm->componentNames[comp]);
1229: PetscStrallocpy(name, &fvm->componentNames[comp]);
1230: return 0;
1231: }
1233: /*@C
1234: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1236: Logically collective on fvm
1237: Input Parameters:
1238: + fvm - the `PetscFV` object
1239: - comp - the component number
1241: Output Parameter:
1242: . name - the component name
1244: Level: intermediate
1246: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1247: @*/
1248: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1249: {
1250: *name = fvm->componentNames[comp];
1251: return 0;
1252: }
1254: /*@
1255: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1257: Logically collective on fvm
1259: Input Parameters:
1260: + fvm - the `PetscFV` object
1261: - dim - The spatial dimension
1263: Level: intermediate
1265: .seealso: `PetscFV`, ``PetscFVGetSpatialDimension()`
1266: @*/
1267: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1268: {
1270: fvm->dim = dim;
1271: return 0;
1272: }
1274: /*@
1275: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1277: Logically collective on fvm
1279: Input Parameter:
1280: . fvm - the `PetscFV` object
1282: Output Parameter:
1283: . dim - The spatial dimension
1285: Level: intermediate
1287: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1288: @*/
1289: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1290: {
1293: *dim = fvm->dim;
1294: return 0;
1295: }
1297: /*@
1298: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1300: Logically collective on fvm
1302: Input Parameters:
1303: + fvm - the `PetscFV` object
1304: - computeGradients - Flag to compute cell gradients
1306: Level: intermediate
1308: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1309: @*/
1310: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1311: {
1313: fvm->computeGradients = computeGradients;
1314: return 0;
1315: }
1317: /*@
1318: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1320: Not collective
1322: Input Parameter:
1323: . fvm - the `PetscFV` object
1325: Output Parameter:
1326: . computeGradients - Flag to compute cell gradients
1328: Level: intermediate
1330: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1331: @*/
1332: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1333: {
1336: *computeGradients = fvm->computeGradients;
1337: return 0;
1338: }
1340: /*@
1341: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1343: Logically collective on fvm
1345: Input Parameters:
1346: + fvm - the `PetscFV` object
1347: - q - The `PetscQuadrature`
1349: Level: intermediate
1351: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1352: @*/
1353: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1354: {
1356: PetscQuadratureDestroy(&fvm->quadrature);
1357: PetscObjectReference((PetscObject)q);
1358: fvm->quadrature = q;
1359: return 0;
1360: }
1362: /*@
1363: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1365: Not collective
1367: Input Parameter:
1368: . fvm - the `PetscFV` object
1370: Output Parameter:
1371: . lim - The `PetscQuadrature`
1373: Level: intermediate
1375: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1376: @*/
1377: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1378: {
1381: if (!fvm->quadrature) {
1382: /* Create default 1-point quadrature */
1383: PetscReal *points, *weights;
1385: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1386: PetscCalloc1(fvm->dim, &points);
1387: PetscMalloc1(1, &weights);
1388: weights[0] = 1.0;
1389: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1390: }
1391: *q = fvm->quadrature;
1392: return 0;
1393: }
1395: /*@
1396: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1398: Not collective
1400: Input Parameter:
1401: . fvm - The `PetscFV` object
1403: Output Parameter:
1404: . sp - The PetscDualSpace object
1406: Level: intermediate
1408: Developer Note:
1409: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1411: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1412: @*/
1413: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1414: {
1417: if (!fvm->dualSpace) {
1418: DM K;
1419: PetscInt dim, Nc, c;
1421: PetscFVGetSpatialDimension(fvm, &dim);
1422: PetscFVGetNumComponents(fvm, &Nc);
1423: PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace);
1424: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1425: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K);
1426: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1427: PetscDualSpaceSetDM(fvm->dualSpace, K);
1428: DMDestroy(&K);
1429: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1430: /* Should we be using PetscFVGetQuadrature() here? */
1431: for (c = 0; c < Nc; ++c) {
1432: PetscQuadrature qc;
1433: PetscReal *points, *weights;
1435: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1436: PetscCalloc1(dim, &points);
1437: PetscCalloc1(Nc, &weights);
1438: weights[c] = 1.0;
1439: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1440: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1441: PetscQuadratureDestroy(&qc);
1442: }
1443: PetscDualSpaceSetUp(fvm->dualSpace);
1444: }
1445: *sp = fvm->dualSpace;
1446: return 0;
1447: }
1449: /*@
1450: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1452: Not collective
1454: Input Parameters:
1455: + fvm - The `PetscFV` object
1456: - sp - The `PetscDualSpace` object
1458: Level: intermediate
1460: Note:
1461: A simple dual space is provided automatically, and the user typically will not need to override it.
1463: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1464: @*/
1465: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1466: {
1469: PetscDualSpaceDestroy(&fvm->dualSpace);
1470: fvm->dualSpace = sp;
1471: PetscObjectReference((PetscObject)fvm->dualSpace);
1472: return 0;
1473: }
1475: /*@C
1476: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1478: Not collective
1480: Input Parameter:
1481: . fvm - The `PetscFV` object
1483: Output Parameter:
1484: . T - The basis function values and derivatives at quadrature points
1486: Level: intermediate
1488: Note:
1489: .vb
1490: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1491: 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
1492: 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
1493: .ve
1495: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1496: @*/
1497: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1498: {
1499: PetscInt npoints;
1500: const PetscReal *points;
1504: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1505: if (!fvm->T) PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);
1506: *T = fvm->T;
1507: return 0;
1508: }
1510: /*@C
1511: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1513: Not collective
1515: Input Parameters:
1516: + fvm - The `PetscFV` object
1517: . nrepl - The number of replicas
1518: . npoints - The number of tabulation points in a replica
1519: . points - The tabulation point coordinates
1520: - K - The order of derivative to tabulate
1522: Output Parameter:
1523: . T - The basis function values and derivative at tabulation points
1525: Level: intermediate
1527: Note:
1528: .vb
1529: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1530: 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
1531: 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
1532: .ve
1534: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1535: @*/
1536: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1537: {
1538: PetscInt pdim = 1; /* Dimension of approximation space P */
1539: PetscInt cdim; /* Spatial dimension */
1540: PetscInt Nc; /* Field components */
1541: PetscInt k, p, d, c, e;
1543: if (!npoints || K < 0) {
1544: *T = NULL;
1545: return 0;
1546: }
1550: PetscFVGetSpatialDimension(fvm, &cdim);
1551: PetscFVGetNumComponents(fvm, &Nc);
1552: PetscMalloc1(1, T);
1553: (*T)->K = !cdim ? 0 : K;
1554: (*T)->Nr = nrepl;
1555: (*T)->Np = npoints;
1556: (*T)->Nb = pdim;
1557: (*T)->Nc = Nc;
1558: (*T)->cdim = cdim;
1559: PetscMalloc1((*T)->K + 1, &(*T)->T);
1560: for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
1561: if (K >= 0) {
1562: for (p = 0; p < nrepl * npoints; ++p)
1563: for (d = 0; d < pdim; ++d)
1564: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1565: }
1566: if (K >= 1) {
1567: for (p = 0; p < nrepl * npoints; ++p)
1568: for (d = 0; d < pdim; ++d)
1569: for (c = 0; c < Nc; ++c)
1570: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1571: }
1572: if (K >= 2) {
1573: for (p = 0; p < nrepl * npoints; ++p)
1574: for (d = 0; d < pdim; ++d)
1575: for (c = 0; c < Nc; ++c)
1576: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1577: }
1578: return 0;
1579: }
1581: /*@C
1582: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1584: Input Parameters:
1585: + fvm - The `PetscFV` object
1586: . numFaces - The number of cell faces which are not constrained
1587: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1589: Level: advanced
1591: .seealso: `PetscFV`, `PetscFVCreate()`
1592: @*/
1593: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1594: {
1596: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1597: return 0;
1598: }
1600: /*@C
1601: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1603: Not collective
1605: Input Parameters:
1606: + fvm - The `PetscFV` object for the field being integrated
1607: . prob - The `PetscDS` specifying the discretizations and continuum functions
1608: . field - The field being integrated
1609: . Nf - The number of faces in the chunk
1610: . fgeom - The face geometry for each face in the chunk
1611: . neighborVol - The volume for each pair of cells in the chunk
1612: . uL - The state from the cell on the left
1613: - uR - The state from the cell on the right
1615: Output Parameters:
1616: + fluxL - the left fluxes for each face
1617: - fluxR - the right fluxes for each face
1619: Level: developer
1621: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1622: @*/
1623: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1624: {
1626: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1627: return 0;
1628: }
1630: /*@
1631: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into smaller copies. This is typically used
1632: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1633: sparsity). It is also used to create an interpolation between regularly refined meshes.
1635: Input Parameter:
1636: . fv - The initial `PetscFV`
1638: Output Parameter:
1639: . fvRef - The refined `PetscFV`
1641: Level: advanced
1643: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1644: @*/
1645: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1646: {
1647: PetscDualSpace Q, Qref;
1648: DM K, Kref;
1649: PetscQuadrature q, qref;
1650: DMPolytopeType ct;
1651: DMPlexTransform tr;
1652: PetscReal *v0;
1653: PetscReal *jac, *invjac;
1654: PetscInt numComp, numSubelements, s;
1656: PetscFVGetDualSpace(fv, &Q);
1657: PetscFVGetQuadrature(fv, &q);
1658: PetscDualSpaceGetDM(Q, &K);
1659: /* Create dual space */
1660: PetscDualSpaceDuplicate(Q, &Qref);
1661: DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref);
1662: PetscDualSpaceSetDM(Qref, Kref);
1663: DMDestroy(&Kref);
1664: PetscDualSpaceSetUp(Qref);
1665: /* Create volume */
1666: PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef);
1667: PetscFVSetDualSpace(*fvRef, Qref);
1668: PetscFVGetNumComponents(fv, &numComp);
1669: PetscFVSetNumComponents(*fvRef, numComp);
1670: PetscFVSetUp(*fvRef);
1671: /* Create quadrature */
1672: DMPlexGetCellType(K, 0, &ct);
1673: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1674: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1675: DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1676: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1677: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1678: for (s = 0; s < numSubelements; ++s) {
1679: PetscQuadrature qs;
1680: const PetscReal *points, *weights;
1681: PetscReal *p, *w;
1682: PetscInt dim, Nc, npoints, np;
1684: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1685: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1686: np = npoints / numSubelements;
1687: PetscMalloc1(np * dim, &p);
1688: PetscMalloc1(np * Nc, &w);
1689: PetscArraycpy(p, &points[s * np * dim], np * dim);
1690: PetscArraycpy(w, &weights[s * np * Nc], np * Nc);
1691: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1692: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1693: PetscQuadratureDestroy(&qs);
1694: }
1695: PetscFVSetQuadrature(*fvRef, qref);
1696: DMPlexTransformDestroy(&tr);
1697: PetscQuadratureDestroy(&qref);
1698: PetscDualSpaceDestroy(&Qref);
1699: return 0;
1700: }
1702: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1703: {
1704: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1706: PetscFree(b);
1707: return 0;
1708: }
1710: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1711: {
1712: PetscInt Nc, c;
1713: PetscViewerFormat format;
1715: PetscFVGetNumComponents(fv, &Nc);
1716: PetscViewerGetFormat(viewer, &format);
1717: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1718: PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc);
1719: for (c = 0; c < Nc; c++) {
1720: if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1721: }
1722: return 0;
1723: }
1725: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1726: {
1727: PetscBool iascii;
1731: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1732: if (iascii) PetscFVView_Upwind_Ascii(fv, viewer);
1733: return 0;
1734: }
1736: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1737: {
1738: return 0;
1739: }
1741: /*
1742: neighborVol[f*2+0] contains the left geom
1743: neighborVol[f*2+1] contains the right geom
1744: */
1745: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1746: {
1747: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1748: void *rctx;
1749: PetscScalar *flux = fvm->fluxWork;
1750: const PetscScalar *constants;
1751: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1753: PetscDSGetTotalComponents(prob, &Nc);
1754: PetscDSGetTotalDimension(prob, &totDim);
1755: PetscDSGetFieldOffset(prob, field, &off);
1756: PetscDSGetRiemannSolver(prob, field, &riemann);
1757: PetscDSGetContext(prob, field, &rctx);
1758: PetscDSGetConstants(prob, &numConstants, &constants);
1759: PetscFVGetSpatialDimension(fvm, &dim);
1760: PetscFVGetNumComponents(fvm, &pdim);
1761: for (f = 0; f < Nf; ++f) {
1762: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1763: for (d = 0; d < pdim; ++d) {
1764: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1765: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1766: }
1767: }
1768: return 0;
1769: }
1771: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1772: {
1773: fvm->ops->setfromoptions = NULL;
1774: fvm->ops->setup = PetscFVSetUp_Upwind;
1775: fvm->ops->view = PetscFVView_Upwind;
1776: fvm->ops->destroy = PetscFVDestroy_Upwind;
1777: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1778: return 0;
1779: }
1781: /*MC
1782: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1784: Level: intermediate
1786: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1787: M*/
1789: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1790: {
1791: PetscFV_Upwind *b;
1794: PetscNew(&b);
1795: fvm->data = b;
1797: PetscFVInitialize_Upwind(fvm);
1798: return 0;
1799: }
1801: #include <petscblaslapack.h>
1803: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1804: {
1805: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1807: PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1808: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1809: PetscFree(ls);
1810: return 0;
1811: }
1813: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1814: {
1815: PetscInt Nc, c;
1816: PetscViewerFormat format;
1818: PetscFVGetNumComponents(fv, &Nc);
1819: PetscViewerGetFormat(viewer, &format);
1820: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1821: PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc);
1822: for (c = 0; c < Nc; c++) {
1823: if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1824: }
1825: return 0;
1826: }
1828: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1829: {
1830: PetscBool iascii;
1834: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1835: if (iascii) PetscFVView_LeastSquares_Ascii(fv, viewer);
1836: return 0;
1837: }
1839: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1840: {
1841: return 0;
1842: }
1844: /* Overwrites A. Can only handle full-rank problems with m>=n */
1845: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1846: {
1847: PetscBool debug = PETSC_FALSE;
1848: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
1849: PetscScalar *R, *Q, *Aback, Alpha;
1851: if (debug) {
1852: PetscMalloc1(m * n, &Aback);
1853: PetscArraycpy(Aback, A, m * n);
1854: }
1856: PetscBLASIntCast(m, &M);
1857: PetscBLASIntCast(n, &N);
1858: PetscBLASIntCast(mstride, &lda);
1859: PetscBLASIntCast(worksize, &ldwork);
1860: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1861: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1862: PetscFPTrapPop();
1864: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1866: /* Extract an explicit representation of Q */
1867: Q = Ainv;
1868: PetscArraycpy(Q, A, mstride * n);
1869: K = N; /* full rank */
1870: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
1873: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1874: Alpha = 1.0;
1875: ldb = lda;
1876: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1877: /* Ainv is Q, overwritten with inverse */
1879: if (debug) { /* Check that pseudo-inverse worked */
1880: PetscScalar Beta = 0.0;
1881: PetscBLASInt ldc;
1882: K = N;
1883: ldc = N;
1884: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1885: PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF);
1886: PetscFree(Aback);
1887: }
1888: return 0;
1889: }
1891: /* Overwrites A. Can handle degenerate problems and m<n. */
1892: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1893: {
1894: PetscScalar *Brhs;
1895: PetscScalar *tmpwork;
1896: PetscReal rcond;
1897: #if defined(PETSC_USE_COMPLEX)
1898: PetscInt rworkSize;
1899: PetscReal *rwork;
1900: #endif
1901: PetscInt i, j, maxmn;
1902: PetscBLASInt M, N, lda, ldb, ldwork;
1903: PetscBLASInt nrhs, irank, info;
1905: /* initialize to identity */
1906: tmpwork = work;
1907: Brhs = Ainv;
1908: maxmn = PetscMax(m, n);
1909: for (j = 0; j < maxmn; j++) {
1910: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
1911: }
1913: PetscBLASIntCast(m, &M);
1914: PetscBLASIntCast(n, &N);
1915: PetscBLASIntCast(mstride, &lda);
1916: PetscBLASIntCast(maxmn, &ldb);
1917: PetscBLASIntCast(worksize, &ldwork);
1918: rcond = -1;
1919: nrhs = M;
1920: #if defined(PETSC_USE_COMPLEX)
1921: rworkSize = 5 * PetscMin(M, N);
1922: PetscMalloc1(rworkSize, &rwork);
1923: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1924: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
1925: PetscFPTrapPop();
1926: PetscFree(rwork);
1927: #else
1928: nrhs = M;
1929: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1930: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
1931: PetscFPTrapPop();
1932: #endif
1934: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1936: return 0;
1937: }
1939: #if 0
1940: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1941: {
1942: PetscReal grad[2] = {0, 0};
1943: const PetscInt *faces;
1944: PetscInt numFaces, f;
1946: DMPlexGetConeSize(dm, cell, &numFaces);
1947: DMPlexGetCone(dm, cell, &faces);
1948: for (f = 0; f < numFaces; ++f) {
1949: const PetscInt *fcells;
1950: const CellGeom *cg1;
1951: const FaceGeom *fg;
1953: DMPlexGetSupport(dm, faces[f], &fcells);
1954: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1955: for (i = 0; i < 2; ++i) {
1956: PetscScalar du;
1958: if (fcells[i] == c) continue;
1959: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1960: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1961: grad[0] += fg->grad[!i][0] * du;
1962: grad[1] += fg->grad[!i][1] * du;
1963: }
1964: }
1965: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1966: return 0;
1967: }
1968: #endif
1970: /*
1971: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1973: Input Parameters:
1974: + fvm - The `PetscFV` object
1975: . numFaces - The number of cell faces which are not constrained
1976: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1978: Level: developer
1980: .seealso: `PetscFV`, `PetscFVCreate()`
1981: */
1982: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1983: {
1984: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1985: const PetscBool useSVD = PETSC_TRUE;
1986: const PetscInt maxFaces = ls->maxFaces;
1987: PetscInt dim, f, d;
1989: if (numFaces > maxFaces) {
1991: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
1992: }
1993: PetscFVGetSpatialDimension(fvm, &dim);
1994: for (f = 0; f < numFaces; ++f) {
1995: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
1996: }
1997: /* Overwrites B with garbage, returns Binv in row-major format */
1998: if (useSVD) {
1999: PetscInt maxmn = PetscMax(numFaces, dim);
2000: PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2001: /* Binv shaped in column-major, coldim=maxmn.*/
2002: for (f = 0; f < numFaces; ++f) {
2003: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2004: }
2005: } else {
2006: PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2007: /* Binv shaped in row-major, rowdim=maxFaces.*/
2008: for (f = 0; f < numFaces; ++f) {
2009: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2010: }
2011: }
2012: return 0;
2013: }
2015: /*
2016: neighborVol[f*2+0] contains the left geom
2017: neighborVol[f*2+1] contains the right geom
2018: */
2019: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2020: {
2021: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2022: void *rctx;
2023: PetscScalar *flux = fvm->fluxWork;
2024: const PetscScalar *constants;
2025: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2027: PetscDSGetTotalComponents(prob, &Nc);
2028: PetscDSGetTotalDimension(prob, &totDim);
2029: PetscDSGetFieldOffset(prob, field, &off);
2030: PetscDSGetRiemannSolver(prob, field, &riemann);
2031: PetscDSGetContext(prob, field, &rctx);
2032: PetscDSGetConstants(prob, &numConstants, &constants);
2033: PetscFVGetSpatialDimension(fvm, &dim);
2034: PetscFVGetNumComponents(fvm, &pdim);
2035: for (f = 0; f < Nf; ++f) {
2036: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2037: for (d = 0; d < pdim; ++d) {
2038: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2039: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2040: }
2041: }
2042: return 0;
2043: }
2045: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2046: {
2047: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2048: PetscInt dim, m, n, nrhs, minmn, maxmn;
2051: PetscFVGetSpatialDimension(fvm, &dim);
2052: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2053: ls->maxFaces = maxFaces;
2054: m = ls->maxFaces;
2055: n = dim;
2056: nrhs = ls->maxFaces;
2057: minmn = PetscMin(m, n);
2058: maxmn = PetscMax(m, n);
2059: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2060: PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work);
2061: return 0;
2062: }
2064: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2065: {
2066: fvm->ops->setfromoptions = NULL;
2067: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2068: fvm->ops->view = PetscFVView_LeastSquares;
2069: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2070: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2071: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2072: return 0;
2073: }
2075: /*MC
2076: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2078: Level: intermediate
2080: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2081: M*/
2083: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2084: {
2085: PetscFV_LeastSquares *ls;
2088: PetscNew(&ls);
2089: fvm->data = ls;
2091: ls->maxFaces = -1;
2092: ls->workSize = -1;
2093: ls->B = NULL;
2094: ls->Binv = NULL;
2095: ls->tau = NULL;
2096: ls->work = NULL;
2098: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2099: PetscFVInitialize_LeastSquares(fvm);
2100: PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2101: return 0;
2102: }
2104: /*@
2105: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2107: Not collective
2109: Input parameters:
2110: + fvm - The `PetscFV` object
2111: - maxFaces - The maximum number of cell faces
2113: Level: intermediate
2115: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2116: @*/
2117: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2118: {
2120: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2121: return 0;
2122: }