Actual source code: fv.c
1: #include <petsc/private/petscfvimpl.h>
2: #include <petscdmplex.h>
3: #include <petscdmplextransform.h>
4: #include <petscds.h>
6: PetscClassId PETSCLIMITER_CLASSID = 0;
8: PetscFunctionList PetscLimiterList = NULL;
9: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
11: PetscBool Limitercite = PETSC_FALSE;
12: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
13: " title = {Analysis of slope limiters on irregular grids},\n"
14: " journal = {AIAA paper},\n"
15: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
16: " volume = {490},\n"
17: " year = {2005}\n}\n";
19: /*@C
20: PetscLimiterRegister - Adds a new `PetscLimiter` implementation
22: Not Collective
24: Input Parameters:
25: + sname - The name of a new user-defined creation routine
26: - function - The creation routine
28: Example Usage:
29: .vb
30: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
31: .ve
33: Then, your `PetscLimiter` type can be chosen with the procedural interface via
34: .vb
35: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
36: PetscLimiterSetType(PetscLimiter, "my_lim");
37: .ve
38: or at runtime via the option
39: .vb
40: -petsclimiter_type my_lim
41: .ve
43: Level: advanced
45: Note:
46: `PetscLimiterRegister()` may be called multiple times to add several user-defined PetscLimiters
48: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterRegisterAll()`, `PetscLimiterRegisterDestroy()`
49: @*/
50: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
51: {
52: PetscFunctionBegin;
53: PetscCall(PetscFunctionListAdd(&PetscLimiterList, sname, function));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: /*@C
58: PetscLimiterSetType - Builds a `PetscLimiter` for a given `PetscLimiterType`
60: Collective
62: Input Parameters:
63: + lim - The `PetscLimiter` object
64: - name - The kind of limiter
66: Options Database Key:
67: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
69: Level: intermediate
71: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
72: @*/
73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
74: {
75: PetscErrorCode (*r)(PetscLimiter);
76: PetscBool match;
78: PetscFunctionBegin;
80: PetscCall(PetscObjectTypeCompare((PetscObject)lim, name, &match));
81: if (match) PetscFunctionReturn(PETSC_SUCCESS);
83: PetscCall(PetscLimiterRegisterAll());
84: PetscCall(PetscFunctionListFind(PetscLimiterList, name, &r));
85: PetscCheck(r, PetscObjectComm((PetscObject)lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
87: PetscTryTypeMethod(lim, destroy);
88: lim->ops->destroy = NULL;
90: PetscCall((*r)(lim));
91: PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@C
96: PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.
98: Not Collective
100: Input Parameter:
101: . lim - The `PetscLimiter`
103: Output Parameter:
104: . name - The `PetscLimiterType`
106: Level: intermediate
108: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
112: PetscFunctionBegin;
114: PetscAssertPointer(name, 2);
115: PetscCall(PetscLimiterRegisterAll());
116: *name = ((PetscObject)lim)->type_name;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@C
121: PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database
123: Collective
125: Input Parameters:
126: + A - the `PetscLimiter` object to view
127: . obj - Optional object that provides the options prefix to use
128: - name - command line option name
130: Level: intermediate
132: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
133: @*/
134: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
135: {
136: PetscFunctionBegin;
138: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*@C
143: PetscLimiterView - Views a `PetscLimiter`
145: Collective
147: Input Parameters:
148: + lim - the `PetscLimiter` object to view
149: - v - the viewer
151: Level: beginner
153: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
154: @*/
155: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
156: {
157: PetscFunctionBegin;
159: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
160: PetscTryTypeMethod(lim, view, v);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
167: Collective
169: Input Parameter:
170: . lim - the `PetscLimiter` object to set options for
172: Level: intermediate
174: .seealso: `PetscLimiter`, `PetscLimiterView()`
175: @*/
176: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
177: {
178: const char *defaultType;
179: char name[256];
180: PetscBool flg;
182: PetscFunctionBegin;
184: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
185: else defaultType = ((PetscObject)lim)->type_name;
186: PetscCall(PetscLimiterRegisterAll());
188: PetscObjectOptionsBegin((PetscObject)lim);
189: PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg));
190: if (flg) {
191: PetscCall(PetscLimiterSetType(lim, name));
192: } else if (!((PetscObject)lim)->type_name) {
193: PetscCall(PetscLimiterSetType(lim, defaultType));
194: }
195: PetscTryTypeMethod(lim, setfromoptions);
196: /* process any options handlers added with PetscObjectAddOptionsHandler() */
197: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
198: PetscOptionsEnd();
199: PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@C
204: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
206: Collective
208: Input Parameter:
209: . lim - the `PetscLimiter` object to setup
211: Level: intermediate
213: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscLimiterDestroy()`
214: @*/
215: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
216: {
217: PetscFunctionBegin;
219: PetscTryTypeMethod(lim, setup);
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: PetscLimiterDestroy - Destroys a `PetscLimiter` object
226: Collective
228: Input Parameter:
229: . lim - the `PetscLimiter` object to destroy
231: Level: beginner
233: .seealso: `PetscLimiter`, `PetscLimiterView()`
234: @*/
235: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
236: {
237: PetscFunctionBegin;
238: if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);
241: if (--((PetscObject)*lim)->refct > 0) {
242: *lim = NULL;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
245: ((PetscObject)*lim)->refct = 0;
247: PetscTryTypeMethod(*lim, destroy);
248: PetscCall(PetscHeaderDestroy(lim));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@
253: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
255: Collective
257: Input Parameter:
258: . comm - The communicator for the `PetscLimiter` object
260: Output Parameter:
261: . lim - The `PetscLimiter` object
263: Level: beginner
265: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
266: @*/
267: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
268: {
269: PetscLimiter l;
271: PetscFunctionBegin;
272: PetscAssertPointer(lim, 2);
273: PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
274: *lim = NULL;
275: PetscCall(PetscFVInitializePackage());
277: PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));
279: *lim = l;
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*@
284: PetscLimiterLimit - Limit the flux
286: Input Parameters:
287: + lim - The `PetscLimiter`
288: - flim - The input field
290: Output Parameter:
291: . phi - The limited field
293: Level: beginner
295: Note:
296: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
297: .vb
298: The classical flux-limited formulation is psi(r) where
300: r = (u[0] - u[-1]) / (u[1] - u[0])
302: The second order TVD region is bounded by
304: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
306: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
307: phi(r)(r+1)/2 in which the reconstructed interface values are
309: u(v) = u[0] + phi(r) (grad u)[0] v
311: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
313: 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))
315: For a nicer symmetric formulation, rewrite in terms of
317: f = (u[0] - u[-1]) / (u[1] - u[-1])
319: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
321: phi(r) = phi(1/r)
323: becomes
325: w(f) = w(1-f).
327: The limiters below implement this final form w(f). The reference methods are
329: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
330: .ve
332: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
333: @*/
334: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
335: {
336: PetscFunctionBegin;
338: PetscAssertPointer(phi, 3);
339: PetscUseTypeMethod(lim, limit, flim, phi);
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
344: {
345: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
347: PetscFunctionBegin;
348: PetscCall(PetscFree(l));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
353: {
354: PetscViewerFormat format;
356: PetscFunctionBegin;
357: PetscCall(PetscViewerGetFormat(viewer, &format));
358: PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
363: {
364: PetscBool iascii;
366: PetscFunctionBegin;
369: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
370: if (iascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
375: {
376: PetscFunctionBegin;
377: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
382: {
383: PetscFunctionBegin;
384: lim->ops->view = PetscLimiterView_Sin;
385: lim->ops->destroy = PetscLimiterDestroy_Sin;
386: lim->ops->limit = PetscLimiterLimit_Sin;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*MC
391: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
393: Level: intermediate
395: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
396: M*/
398: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
399: {
400: PetscLimiter_Sin *l;
402: PetscFunctionBegin;
404: PetscCall(PetscNew(&l));
405: lim->data = l;
407: PetscCall(PetscLimiterInitialize_Sin(lim));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
412: {
413: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
415: PetscFunctionBegin;
416: PetscCall(PetscFree(l));
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
421: {
422: PetscViewerFormat format;
424: PetscFunctionBegin;
425: PetscCall(PetscViewerGetFormat(viewer, &format));
426: PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
431: {
432: PetscBool iascii;
434: PetscFunctionBegin;
437: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
438: if (iascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
443: {
444: PetscFunctionBegin;
445: *phi = 0.0;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
450: {
451: PetscFunctionBegin;
452: lim->ops->view = PetscLimiterView_Zero;
453: lim->ops->destroy = PetscLimiterDestroy_Zero;
454: lim->ops->limit = PetscLimiterLimit_Zero;
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: /*MC
459: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
461: Level: intermediate
463: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
464: M*/
466: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
467: {
468: PetscLimiter_Zero *l;
470: PetscFunctionBegin;
472: PetscCall(PetscNew(&l));
473: lim->data = l;
475: PetscCall(PetscLimiterInitialize_Zero(lim));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
480: {
481: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
483: PetscFunctionBegin;
484: PetscCall(PetscFree(l));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
489: {
490: PetscViewerFormat format;
492: PetscFunctionBegin;
493: PetscCall(PetscViewerGetFormat(viewer, &format));
494: PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
499: {
500: PetscBool iascii;
502: PetscFunctionBegin;
505: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
506: if (iascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
511: {
512: PetscFunctionBegin;
513: *phi = 1.0;
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
518: {
519: PetscFunctionBegin;
520: lim->ops->view = PetscLimiterView_None;
521: lim->ops->destroy = PetscLimiterDestroy_None;
522: lim->ops->limit = PetscLimiterLimit_None;
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*MC
527: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
529: Level: intermediate
531: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
532: M*/
534: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
535: {
536: PetscLimiter_None *l;
538: PetscFunctionBegin;
540: PetscCall(PetscNew(&l));
541: lim->data = l;
543: PetscCall(PetscLimiterInitialize_None(lim));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
548: {
549: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
551: PetscFunctionBegin;
552: PetscCall(PetscFree(l));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
557: {
558: PetscViewerFormat format;
560: PetscFunctionBegin;
561: PetscCall(PetscViewerGetFormat(viewer, &format));
562: PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
567: {
568: PetscBool iascii;
570: PetscFunctionBegin;
573: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
574: if (iascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
579: {
580: PetscFunctionBegin;
581: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
586: {
587: PetscFunctionBegin;
588: lim->ops->view = PetscLimiterView_Minmod;
589: lim->ops->destroy = PetscLimiterDestroy_Minmod;
590: lim->ops->limit = PetscLimiterLimit_Minmod;
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: /*MC
595: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
597: Level: intermediate
599: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
600: M*/
602: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
603: {
604: PetscLimiter_Minmod *l;
606: PetscFunctionBegin;
608: PetscCall(PetscNew(&l));
609: lim->data = l;
611: PetscCall(PetscLimiterInitialize_Minmod(lim));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
616: {
617: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
619: PetscFunctionBegin;
620: PetscCall(PetscFree(l));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
625: {
626: PetscViewerFormat format;
628: PetscFunctionBegin;
629: PetscCall(PetscViewerGetFormat(viewer, &format));
630: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
635: {
636: PetscBool iascii;
638: PetscFunctionBegin;
641: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
642: if (iascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
647: {
648: PetscFunctionBegin;
649: *phi = PetscMax(0, 4 * f * (1 - f));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
654: {
655: PetscFunctionBegin;
656: lim->ops->view = PetscLimiterView_VanLeer;
657: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
658: lim->ops->limit = PetscLimiterLimit_VanLeer;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*MC
663: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
665: Level: intermediate
667: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
668: M*/
670: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
671: {
672: PetscLimiter_VanLeer *l;
674: PetscFunctionBegin;
676: PetscCall(PetscNew(&l));
677: lim->data = l;
679: PetscCall(PetscLimiterInitialize_VanLeer(lim));
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
684: {
685: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
687: PetscFunctionBegin;
688: PetscCall(PetscFree(l));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
693: {
694: PetscViewerFormat format;
696: PetscFunctionBegin;
697: PetscCall(PetscViewerGetFormat(viewer, &format));
698: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
703: {
704: PetscBool iascii;
706: PetscFunctionBegin;
709: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
710: if (iascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
715: {
716: PetscFunctionBegin;
717: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
722: {
723: PetscFunctionBegin;
724: lim->ops->view = PetscLimiterView_VanAlbada;
725: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
726: lim->ops->limit = PetscLimiterLimit_VanAlbada;
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: /*MC
731: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
733: Level: intermediate
735: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
736: M*/
738: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
739: {
740: PetscLimiter_VanAlbada *l;
742: PetscFunctionBegin;
744: PetscCall(PetscNew(&l));
745: lim->data = l;
747: PetscCall(PetscLimiterInitialize_VanAlbada(lim));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
752: {
753: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
755: PetscFunctionBegin;
756: PetscCall(PetscFree(l));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
761: {
762: PetscViewerFormat format;
764: PetscFunctionBegin;
765: PetscCall(PetscViewerGetFormat(viewer, &format));
766: PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
767: PetscFunctionReturn(PETSC_SUCCESS);
768: }
770: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
771: {
772: PetscBool iascii;
774: PetscFunctionBegin;
777: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
778: if (iascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
783: {
784: PetscFunctionBegin;
785: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
790: {
791: PetscFunctionBegin;
792: lim->ops->view = PetscLimiterView_Superbee;
793: lim->ops->destroy = PetscLimiterDestroy_Superbee;
794: lim->ops->limit = PetscLimiterLimit_Superbee;
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*MC
799: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
801: Level: intermediate
803: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
804: M*/
806: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
807: {
808: PetscLimiter_Superbee *l;
810: PetscFunctionBegin;
812: PetscCall(PetscNew(&l));
813: lim->data = l;
815: PetscCall(PetscLimiterInitialize_Superbee(lim));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
820: {
821: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
823: PetscFunctionBegin;
824: PetscCall(PetscFree(l));
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }
828: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
829: {
830: PetscViewerFormat format;
832: PetscFunctionBegin;
833: PetscCall(PetscViewerGetFormat(viewer, &format));
834: PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
839: {
840: PetscBool iascii;
842: PetscFunctionBegin;
845: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
846: if (iascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /* aka Barth-Jespersen */
851: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
852: {
853: PetscFunctionBegin;
854: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
859: {
860: PetscFunctionBegin;
861: lim->ops->view = PetscLimiterView_MC;
862: lim->ops->destroy = PetscLimiterDestroy_MC;
863: lim->ops->limit = PetscLimiterLimit_MC;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*MC
868: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
870: Level: intermediate
872: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
873: M*/
875: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
876: {
877: PetscLimiter_MC *l;
879: PetscFunctionBegin;
881: PetscCall(PetscNew(&l));
882: lim->data = l;
884: PetscCall(PetscLimiterInitialize_MC(lim));
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
888: PetscClassId PETSCFV_CLASSID = 0;
890: PetscFunctionList PetscFVList = NULL;
891: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
893: /*@C
894: PetscFVRegister - Adds a new `PetscFV` implementation
896: Not Collective
898: Input Parameters:
899: + sname - The name of a new user-defined creation routine
900: - function - The creation routine itself
902: Example Usage:
903: .vb
904: PetscFVRegister("my_fv", MyPetscFVCreate);
905: .ve
907: Then, your PetscFV type can be chosen with the procedural interface via
908: .vb
909: PetscFVCreate(MPI_Comm, PetscFV *);
910: PetscFVSetType(PetscFV, "my_fv");
911: .ve
912: or at runtime via the option
913: .vb
914: -petscfv_type my_fv
915: .ve
917: Level: advanced
919: Note:
920: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
922: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
923: @*/
924: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
925: {
926: PetscFunctionBegin;
927: PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: /*@C
932: PetscFVSetType - Builds a particular `PetscFV`
934: Collective
936: Input Parameters:
937: + fvm - The `PetscFV` object
938: - name - The type of FVM space
940: Options Database Key:
941: . -petscfv_type <type> - Sets the `PetscFVType`; use -help for a list of available types
943: Level: intermediate
945: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
946: @*/
947: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
948: {
949: PetscErrorCode (*r)(PetscFV);
950: PetscBool match;
952: PetscFunctionBegin;
954: PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
955: if (match) PetscFunctionReturn(PETSC_SUCCESS);
957: PetscCall(PetscFVRegisterAll());
958: PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
959: PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
961: PetscTryTypeMethod(fvm, destroy);
962: fvm->ops->destroy = NULL;
964: PetscCall((*r)(fvm));
965: PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: /*@C
970: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
972: Not Collective
974: Input Parameter:
975: . fvm - The `PetscFV`
977: Output Parameter:
978: . name - The `PetscFVType` name
980: Level: intermediate
982: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
983: @*/
984: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
985: {
986: PetscFunctionBegin;
988: PetscAssertPointer(name, 2);
989: PetscCall(PetscFVRegisterAll());
990: *name = ((PetscObject)fvm)->type_name;
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: /*@C
995: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
997: Collective
999: Input Parameters:
1000: + A - the `PetscFV` object
1001: . obj - Optional object that provides the options prefix
1002: - name - command line option name
1004: Level: intermediate
1006: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
1007: @*/
1008: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1009: {
1010: PetscFunctionBegin;
1012: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: /*@C
1017: PetscFVView - Views a `PetscFV`
1019: Collective
1021: Input Parameters:
1022: + fvm - the `PetscFV` object to view
1023: - v - the viewer
1025: Level: beginner
1027: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1028: @*/
1029: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1030: {
1031: PetscFunctionBegin;
1033: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1034: PetscTryTypeMethod(fvm, view, v);
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
1041: Collective
1043: Input Parameter:
1044: . fvm - the `PetscFV` object to set options for
1046: Options Database Key:
1047: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1049: Level: intermediate
1051: .seealso: `PetscFV`, `PetscFVView()`
1052: @*/
1053: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1054: {
1055: const char *defaultType;
1056: char name[256];
1057: PetscBool flg;
1059: PetscFunctionBegin;
1061: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1062: else defaultType = ((PetscObject)fvm)->type_name;
1063: PetscCall(PetscFVRegisterAll());
1065: PetscObjectOptionsBegin((PetscObject)fvm);
1066: PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg));
1067: if (flg) {
1068: PetscCall(PetscFVSetType(fvm, name));
1069: } else if (!((PetscObject)fvm)->type_name) {
1070: PetscCall(PetscFVSetType(fvm, defaultType));
1071: }
1072: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1073: PetscTryTypeMethod(fvm, setfromoptions);
1074: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1075: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1076: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1077: PetscOptionsEnd();
1078: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }
1082: /*@
1083: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1085: Collective
1087: Input Parameter:
1088: . fvm - the `PetscFV` object to setup
1090: Level: intermediate
1092: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1093: @*/
1094: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1095: {
1096: PetscFunctionBegin;
1098: PetscCall(PetscLimiterSetUp(fvm->limiter));
1099: PetscTryTypeMethod(fvm, setup);
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: /*@
1104: PetscFVDestroy - Destroys a `PetscFV` object
1106: Collective
1108: Input Parameter:
1109: . fvm - the `PetscFV` object to destroy
1111: Level: beginner
1113: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1114: @*/
1115: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1116: {
1117: PetscInt i;
1119: PetscFunctionBegin;
1120: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1123: if (--((PetscObject)*fvm)->refct > 0) {
1124: *fvm = NULL;
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1127: ((PetscObject)*fvm)->refct = 0;
1129: for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1130: PetscCall(PetscFree((*fvm)->componentNames));
1131: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1132: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1133: PetscCall(PetscFree((*fvm)->fluxWork));
1134: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1135: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1137: PetscTryTypeMethod(*fvm, destroy);
1138: PetscCall(PetscHeaderDestroy(fvm));
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1142: /*@
1143: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1145: Collective
1147: Input Parameter:
1148: . comm - The communicator for the `PetscFV` object
1150: Output Parameter:
1151: . fvm - The `PetscFV` object
1153: Level: beginner
1155: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1156: @*/
1157: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1158: {
1159: PetscFV f;
1161: PetscFunctionBegin;
1162: PetscAssertPointer(fvm, 2);
1163: *fvm = NULL;
1164: PetscCall(PetscFVInitializePackage());
1166: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1167: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1169: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1170: f->numComponents = 1;
1171: f->dim = 0;
1172: f->computeGradients = PETSC_FALSE;
1173: f->fluxWork = NULL;
1174: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1176: *fvm = f;
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: /*@
1181: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1183: Logically Collective
1185: Input Parameters:
1186: + fvm - the `PetscFV` object
1187: - lim - The `PetscLimiter`
1189: Level: intermediate
1191: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1192: @*/
1193: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1194: {
1195: PetscFunctionBegin;
1198: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1199: PetscCall(PetscObjectReference((PetscObject)lim));
1200: fvm->limiter = lim;
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1207: Not Collective
1209: Input Parameter:
1210: . fvm - the `PetscFV` object
1212: Output Parameter:
1213: . lim - The `PetscLimiter`
1215: Level: intermediate
1217: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1218: @*/
1219: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1220: {
1221: PetscFunctionBegin;
1223: PetscAssertPointer(lim, 2);
1224: *lim = fvm->limiter;
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: /*@
1229: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1231: Logically Collective
1233: Input Parameters:
1234: + fvm - the `PetscFV` object
1235: - comp - The number of components
1237: Level: intermediate
1239: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1240: @*/
1241: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1242: {
1243: PetscFunctionBegin;
1245: if (fvm->numComponents != comp) {
1246: PetscInt i;
1248: for (i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1249: PetscCall(PetscFree(fvm->componentNames));
1250: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1251: }
1252: fvm->numComponents = comp;
1253: PetscCall(PetscFree(fvm->fluxWork));
1254: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: /*@
1259: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1261: Not Collective
1263: Input Parameter:
1264: . fvm - the `PetscFV` object
1266: Output Parameter:
1267: . comp - The number of components
1269: Level: intermediate
1271: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1272: @*/
1273: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1274: {
1275: PetscFunctionBegin;
1277: PetscAssertPointer(comp, 2);
1278: *comp = fvm->numComponents;
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*@C
1283: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1285: Logically Collective
1287: Input Parameters:
1288: + fvm - the `PetscFV` object
1289: . comp - the component number
1290: - name - the component name
1292: Level: intermediate
1294: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1295: @*/
1296: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1297: {
1298: PetscFunctionBegin;
1299: PetscCall(PetscFree(fvm->componentNames[comp]));
1300: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1301: PetscFunctionReturn(PETSC_SUCCESS);
1302: }
1304: /*@C
1305: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1307: 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;
1367: PetscAssertPointer(dim, 2);
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;
1412: PetscAssertPointer(computeGradients, 2);
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(PetscObjectReference((PetscObject)q));
1435: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
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: . q - The `PetscQuadrature`
1451: Level: intermediate
1453: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1454: @*/
1455: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1456: {
1457: PetscFunctionBegin;
1459: PetscAssertPointer(q, 2);
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: PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`
1477: Not Collective
1479: Input Parameters:
1480: + fvm - The `PetscFV` object
1481: - ct - The `DMPolytopeType` for the cell
1483: Level: intermediate
1485: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1486: @*/
1487: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1488: {
1489: DM K;
1490: PetscInt dim, Nc;
1492: PetscFunctionBegin;
1493: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1494: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1495: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1496: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1497: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1498: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1499: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1500: PetscCall(DMDestroy(&K));
1501: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1502: // Should we be using PetscFVGetQuadrature() here?
1503: for (PetscInt c = 0; c < Nc; ++c) {
1504: PetscQuadrature qc;
1505: PetscReal *points, *weights;
1507: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1508: PetscCall(PetscCalloc1(dim, &points));
1509: PetscCall(PetscCalloc1(Nc, &weights));
1510: weights[c] = 1.0;
1511: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1512: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1513: PetscCall(PetscQuadratureDestroy(&qc));
1514: }
1515: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1516: PetscFunctionReturn(PETSC_SUCCESS);
1517: }
1519: /*@
1520: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1522: Not Collective
1524: Input Parameter:
1525: . fvm - The `PetscFV` object
1527: Output Parameter:
1528: . sp - The `PetscDualSpace` object
1530: Level: intermediate
1532: Developer Notes:
1533: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1535: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1536: @*/
1537: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1538: {
1539: PetscFunctionBegin;
1541: PetscAssertPointer(sp, 2);
1542: if (!fvm->dualSpace) {
1543: PetscInt dim;
1545: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1546: PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1547: }
1548: *sp = fvm->dualSpace;
1549: PetscFunctionReturn(PETSC_SUCCESS);
1550: }
1552: /*@
1553: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1555: Not Collective
1557: Input Parameters:
1558: + fvm - The `PetscFV` object
1559: - sp - The `PetscDualSpace` object
1561: Level: intermediate
1563: Note:
1564: A simple dual space is provided automatically, and the user typically will not need to override it.
1566: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1567: @*/
1568: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1569: {
1570: PetscFunctionBegin;
1573: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1574: fvm->dualSpace = sp;
1575: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1576: PetscFunctionReturn(PETSC_SUCCESS);
1577: }
1579: /*@C
1580: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1582: Not Collective
1584: Input Parameter:
1585: . fvm - The `PetscFV` object
1587: Output Parameter:
1588: . T - The basis function values and derivatives at quadrature points
1590: Level: intermediate
1592: Note:
1593: .vb
1594: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1595: 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
1596: 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
1597: .ve
1599: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1600: @*/
1601: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1602: {
1603: PetscInt npoints;
1604: const PetscReal *points;
1606: PetscFunctionBegin;
1608: PetscAssertPointer(T, 2);
1609: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1610: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1611: *T = fvm->T;
1612: PetscFunctionReturn(PETSC_SUCCESS);
1613: }
1615: /*@C
1616: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1618: Not Collective
1620: Input Parameters:
1621: + fvm - The `PetscFV` object
1622: . nrepl - The number of replicas
1623: . npoints - The number of tabulation points in a replica
1624: . points - The tabulation point coordinates
1625: - K - The order of derivative to tabulate
1627: Output Parameter:
1628: . T - The basis function values and derivative at tabulation points
1630: Level: intermediate
1632: Note:
1633: .vb
1634: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1635: 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
1636: 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
1637: .ve
1639: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1640: @*/
1641: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1642: {
1643: PetscInt pdim; // Dimension of approximation space P
1644: PetscInt cdim; // Spatial dimension
1645: PetscInt Nc; // Field components
1646: PetscInt k, p, d, c, e;
1648: PetscFunctionBegin;
1649: if (!npoints || K < 0) {
1650: *T = NULL;
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: PetscAssertPointer(points, 4);
1655: PetscAssertPointer(T, 6);
1656: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1657: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1658: pdim = Nc;
1659: PetscCall(PetscMalloc1(1, T));
1660: (*T)->K = !cdim ? 0 : K;
1661: (*T)->Nr = nrepl;
1662: (*T)->Np = npoints;
1663: (*T)->Nb = pdim;
1664: (*T)->Nc = Nc;
1665: (*T)->cdim = cdim;
1666: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1667: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1668: if (K >= 0) {
1669: for (p = 0; p < nrepl * npoints; ++p)
1670: for (d = 0; d < pdim; ++d)
1671: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1672: }
1673: if (K >= 1) {
1674: for (p = 0; p < nrepl * npoints; ++p)
1675: for (d = 0; d < pdim; ++d)
1676: for (c = 0; c < Nc; ++c)
1677: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1678: }
1679: if (K >= 2) {
1680: for (p = 0; p < nrepl * npoints; ++p)
1681: for (d = 0; d < pdim; ++d)
1682: for (c = 0; c < Nc; ++c)
1683: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1684: }
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@C
1689: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1691: Input Parameters:
1692: + fvm - The `PetscFV` object
1693: . numFaces - The number of cell faces which are not constrained
1694: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1696: Output Parameter:
1697: . grad - the gradient
1699: Level: advanced
1701: .seealso: `PetscFV`, `PetscFVCreate()`
1702: @*/
1703: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1704: {
1705: PetscFunctionBegin;
1707: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: /*@C
1712: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1714: Not Collective
1716: Input Parameters:
1717: + fvm - The `PetscFV` object for the field being integrated
1718: . prob - The `PetscDS` specifying the discretizations and continuum functions
1719: . field - The field being integrated
1720: . Nf - The number of faces in the chunk
1721: . fgeom - The face geometry for each face in the chunk
1722: . neighborVol - The volume for each pair of cells in the chunk
1723: . uL - The state from the cell on the left
1724: - uR - The state from the cell on the right
1726: Output Parameters:
1727: + fluxL - the left fluxes for each face
1728: - fluxR - the right fluxes for each face
1730: Level: developer
1732: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1733: @*/
1734: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1735: {
1736: PetscFunctionBegin;
1738: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.
1745: Input Parameter:
1746: . fv - The initial `PetscFV`
1748: Output Parameter:
1749: . fvNew - A clone of the `PetscFV`
1751: Level: advanced
1753: Notes:
1754: This is typically used to change the number of components.
1756: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1757: @*/
1758: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1759: {
1760: PetscDualSpace Q;
1761: DM K;
1762: PetscQuadrature q;
1763: PetscInt Nc, cdim;
1765: PetscFunctionBegin;
1766: PetscCall(PetscFVGetDualSpace(fv, &Q));
1767: PetscCall(PetscFVGetQuadrature(fv, &q));
1768: PetscCall(PetscDualSpaceGetDM(Q, &K));
1770: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1771: PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1772: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1773: PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1774: PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1775: PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1776: PetscCall(PetscFVSetQuadrature(*fvNew, q));
1777: PetscFunctionReturn(PETSC_SUCCESS);
1778: }
1780: /*@
1781: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1782: smaller copies.
1784: Input Parameter:
1785: . fv - The initial `PetscFV`
1787: Output Parameter:
1788: . fvRef - The refined `PetscFV`
1790: Level: advanced
1792: Notes:
1793: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1794: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1795: interpolation between regularly refined meshes.
1797: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1798: @*/
1799: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1800: {
1801: PetscDualSpace Q, Qref;
1802: DM K, Kref;
1803: PetscQuadrature q, qref;
1804: DMPolytopeType ct;
1805: DMPlexTransform tr;
1806: PetscReal *v0;
1807: PetscReal *jac, *invjac;
1808: PetscInt numComp, numSubelements, s;
1810: PetscFunctionBegin;
1811: PetscCall(PetscFVGetDualSpace(fv, &Q));
1812: PetscCall(PetscFVGetQuadrature(fv, &q));
1813: PetscCall(PetscDualSpaceGetDM(Q, &K));
1814: /* Create dual space */
1815: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1816: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1817: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1818: PetscCall(DMDestroy(&Kref));
1819: PetscCall(PetscDualSpaceSetUp(Qref));
1820: /* Create volume */
1821: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1822: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1823: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1824: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1825: PetscCall(PetscFVSetUp(*fvRef));
1826: /* Create quadrature */
1827: PetscCall(DMPlexGetCellType(K, 0, &ct));
1828: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1829: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1830: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1831: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1832: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1833: for (s = 0; s < numSubelements; ++s) {
1834: PetscQuadrature qs;
1835: const PetscReal *points, *weights;
1836: PetscReal *p, *w;
1837: PetscInt dim, Nc, npoints, np;
1839: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1840: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1841: np = npoints / numSubelements;
1842: PetscCall(PetscMalloc1(np * dim, &p));
1843: PetscCall(PetscMalloc1(np * Nc, &w));
1844: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1845: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1846: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1847: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1848: PetscCall(PetscQuadratureDestroy(&qs));
1849: }
1850: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1851: PetscCall(DMPlexTransformDestroy(&tr));
1852: PetscCall(PetscQuadratureDestroy(&qref));
1853: PetscCall(PetscDualSpaceDestroy(&Qref));
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1858: {
1859: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1861: PetscFunctionBegin;
1862: PetscCall(PetscFree(b));
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1866: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1867: {
1868: PetscInt Nc, c;
1869: PetscViewerFormat format;
1871: PetscFunctionBegin;
1872: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1873: PetscCall(PetscViewerGetFormat(viewer, &format));
1874: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1875: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1876: for (c = 0; c < Nc; c++) {
1877: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1878: }
1879: PetscFunctionReturn(PETSC_SUCCESS);
1880: }
1882: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1883: {
1884: PetscBool iascii;
1886: PetscFunctionBegin;
1889: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1890: if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1891: PetscFunctionReturn(PETSC_SUCCESS);
1892: }
1894: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1895: {
1896: PetscFunctionBegin;
1897: PetscFunctionReturn(PETSC_SUCCESS);
1898: }
1900: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1901: {
1902: PetscInt dim;
1904: PetscFunctionBegin;
1905: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1906: for (PetscInt f = 0; f < numFaces; ++f) {
1907: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1908: }
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*
1913: neighborVol[f*2+0] contains the left geom
1914: neighborVol[f*2+1] contains the right geom
1915: */
1916: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1917: {
1918: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1919: void *rctx;
1920: PetscScalar *flux = fvm->fluxWork;
1921: const PetscScalar *constants;
1922: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1924: PetscFunctionBegin;
1925: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1926: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1927: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1928: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1929: PetscCall(PetscDSGetContext(prob, field, &rctx));
1930: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1931: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1932: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1933: for (f = 0; f < Nf; ++f) {
1934: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1935: for (d = 0; d < pdim; ++d) {
1936: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1937: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1938: }
1939: }
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1944: {
1945: PetscFunctionBegin;
1946: fvm->ops->setfromoptions = NULL;
1947: fvm->ops->setup = PetscFVSetUp_Upwind;
1948: fvm->ops->view = PetscFVView_Upwind;
1949: fvm->ops->destroy = PetscFVDestroy_Upwind;
1950: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1951: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }
1955: /*MC
1956: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1958: Level: intermediate
1960: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1961: M*/
1963: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1964: {
1965: PetscFV_Upwind *b;
1967: PetscFunctionBegin;
1969: PetscCall(PetscNew(&b));
1970: fvm->data = b;
1972: PetscCall(PetscFVInitialize_Upwind(fvm));
1973: PetscFunctionReturn(PETSC_SUCCESS);
1974: }
1976: #include <petscblaslapack.h>
1978: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1979: {
1980: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1982: PetscFunctionBegin;
1983: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1984: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1985: PetscCall(PetscFree(ls));
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1990: {
1991: PetscInt Nc, c;
1992: PetscViewerFormat format;
1994: PetscFunctionBegin;
1995: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1996: PetscCall(PetscViewerGetFormat(viewer, &format));
1997: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1998: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1999: for (c = 0; c < Nc; c++) {
2000: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
2001: }
2002: PetscFunctionReturn(PETSC_SUCCESS);
2003: }
2005: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2006: {
2007: PetscBool iascii;
2009: PetscFunctionBegin;
2012: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2013: if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2014: PetscFunctionReturn(PETSC_SUCCESS);
2015: }
2017: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2018: {
2019: PetscFunctionBegin;
2020: PetscFunctionReturn(PETSC_SUCCESS);
2021: }
2023: /* Overwrites A. Can only handle full-rank problems with m>=n */
2024: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2025: {
2026: PetscBool debug = PETSC_FALSE;
2027: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2028: PetscScalar *R, *Q, *Aback, Alpha;
2030: PetscFunctionBegin;
2031: if (debug) {
2032: PetscCall(PetscMalloc1(m * n, &Aback));
2033: PetscCall(PetscArraycpy(Aback, A, m * n));
2034: }
2036: PetscCall(PetscBLASIntCast(m, &M));
2037: PetscCall(PetscBLASIntCast(n, &N));
2038: PetscCall(PetscBLASIntCast(mstride, &lda));
2039: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2040: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2041: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2042: PetscCall(PetscFPTrapPop());
2043: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2044: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2046: /* Extract an explicit representation of Q */
2047: Q = Ainv;
2048: PetscCall(PetscArraycpy(Q, A, mstride * n));
2049: K = N; /* full rank */
2050: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2051: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2053: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2054: Alpha = 1.0;
2055: ldb = lda;
2056: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2057: /* Ainv is Q, overwritten with inverse */
2059: if (debug) { /* Check that pseudo-inverse worked */
2060: PetscScalar Beta = 0.0;
2061: PetscBLASInt ldc;
2062: K = N;
2063: ldc = N;
2064: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2065: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2066: PetscCall(PetscFree(Aback));
2067: }
2068: PetscFunctionReturn(PETSC_SUCCESS);
2069: }
2071: /* Overwrites A. Can handle degenerate problems and m<n. */
2072: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2073: {
2074: PetscScalar *Brhs;
2075: PetscScalar *tmpwork;
2076: PetscReal rcond;
2077: #if defined(PETSC_USE_COMPLEX)
2078: PetscInt rworkSize;
2079: PetscReal *rwork;
2080: #endif
2081: PetscInt i, j, maxmn;
2082: PetscBLASInt M, N, lda, ldb, ldwork;
2083: PetscBLASInt nrhs, irank, info;
2085: PetscFunctionBegin;
2086: /* initialize to identity */
2087: tmpwork = work;
2088: Brhs = Ainv;
2089: maxmn = PetscMax(m, n);
2090: for (j = 0; j < maxmn; j++) {
2091: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2092: }
2094: PetscCall(PetscBLASIntCast(m, &M));
2095: PetscCall(PetscBLASIntCast(n, &N));
2096: PetscCall(PetscBLASIntCast(mstride, &lda));
2097: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2098: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2099: rcond = -1;
2100: nrhs = M;
2101: #if defined(PETSC_USE_COMPLEX)
2102: rworkSize = 5 * PetscMin(M, N);
2103: PetscCall(PetscMalloc1(rworkSize, &rwork));
2104: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2105: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2106: PetscCall(PetscFPTrapPop());
2107: PetscCall(PetscFree(rwork));
2108: #else
2109: nrhs = M;
2110: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2111: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2112: PetscCall(PetscFPTrapPop());
2113: #endif
2114: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2115: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2116: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }
2120: #if 0
2121: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2122: {
2123: PetscReal grad[2] = {0, 0};
2124: const PetscInt *faces;
2125: PetscInt numFaces, f;
2127: PetscFunctionBegin;
2128: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2129: PetscCall(DMPlexGetCone(dm, cell, &faces));
2130: for (f = 0; f < numFaces; ++f) {
2131: const PetscInt *fcells;
2132: const CellGeom *cg1;
2133: const FaceGeom *fg;
2135: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2136: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2137: for (i = 0; i < 2; ++i) {
2138: PetscScalar du;
2140: if (fcells[i] == c) continue;
2141: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2142: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2143: grad[0] += fg->grad[!i][0] * du;
2144: grad[1] += fg->grad[!i][1] * du;
2145: }
2146: }
2147: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2148: PetscFunctionReturn(PETSC_SUCCESS);
2149: }
2150: #endif
2152: /*
2153: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2155: Input Parameters:
2156: + fvm - The `PetscFV` object
2157: . numFaces - The number of cell faces which are not constrained
2158: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2160: Level: developer
2162: .seealso: `PetscFV`, `PetscFVCreate()`
2163: */
2164: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2165: {
2166: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2167: const PetscBool useSVD = PETSC_TRUE;
2168: const PetscInt maxFaces = ls->maxFaces;
2169: PetscInt dim, f, d;
2171: PetscFunctionBegin;
2172: if (numFaces > maxFaces) {
2173: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2174: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2175: }
2176: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2177: for (f = 0; f < numFaces; ++f) {
2178: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2179: }
2180: /* Overwrites B with garbage, returns Binv in row-major format */
2181: if (useSVD) {
2182: PetscInt maxmn = PetscMax(numFaces, dim);
2183: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2184: /* Binv shaped in column-major, coldim=maxmn.*/
2185: for (f = 0; f < numFaces; ++f) {
2186: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2187: }
2188: } else {
2189: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2190: /* Binv shaped in row-major, rowdim=maxFaces.*/
2191: for (f = 0; f < numFaces; ++f) {
2192: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2193: }
2194: }
2195: PetscFunctionReturn(PETSC_SUCCESS);
2196: }
2198: /*
2199: neighborVol[f*2+0] contains the left geom
2200: neighborVol[f*2+1] contains the right geom
2201: */
2202: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2203: {
2204: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2205: void *rctx;
2206: PetscScalar *flux = fvm->fluxWork;
2207: const PetscScalar *constants;
2208: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2210: PetscFunctionBegin;
2211: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2212: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2213: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2214: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2215: PetscCall(PetscDSGetContext(prob, field, &rctx));
2216: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2217: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2218: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2219: for (f = 0; f < Nf; ++f) {
2220: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2221: for (d = 0; d < pdim; ++d) {
2222: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2223: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2224: }
2225: }
2226: PetscFunctionReturn(PETSC_SUCCESS);
2227: }
2229: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2230: {
2231: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2232: PetscInt dim, m, n, nrhs, minmn, maxmn;
2234: PetscFunctionBegin;
2236: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2237: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2238: ls->maxFaces = maxFaces;
2239: m = ls->maxFaces;
2240: n = dim;
2241: nrhs = ls->maxFaces;
2242: minmn = PetscMin(m, n);
2243: maxmn = PetscMax(m, n);
2244: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2245: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2246: PetscFunctionReturn(PETSC_SUCCESS);
2247: }
2249: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2250: {
2251: PetscFunctionBegin;
2252: fvm->ops->setfromoptions = NULL;
2253: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2254: fvm->ops->view = PetscFVView_LeastSquares;
2255: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2256: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2257: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2258: PetscFunctionReturn(PETSC_SUCCESS);
2259: }
2261: /*MC
2262: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2264: Level: intermediate
2266: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2267: M*/
2269: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2270: {
2271: PetscFV_LeastSquares *ls;
2273: PetscFunctionBegin;
2275: PetscCall(PetscNew(&ls));
2276: fvm->data = ls;
2278: ls->maxFaces = -1;
2279: ls->workSize = -1;
2280: ls->B = NULL;
2281: ls->Binv = NULL;
2282: ls->tau = NULL;
2283: ls->work = NULL;
2285: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2286: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2287: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2288: PetscFunctionReturn(PETSC_SUCCESS);
2289: }
2291: /*@
2292: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2294: Not Collective
2296: Input Parameters:
2297: + fvm - The `PetscFV` object
2298: - maxFaces - The maximum number of cell faces
2300: Level: intermediate
2302: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2303: @*/
2304: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2305: {
2306: PetscFunctionBegin;
2308: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2309: PetscFunctionReturn(PETSC_SUCCESS);
2310: }