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, No Fortran Support
24: Input Parameters:
25: + sname - The name of a new user-defined creation routine
26: - function - The creation routine
28: Example Usage:
29: .vb
30: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
31: .ve
33: Then, your `PetscLimiter` type can be chosen with the procedural interface via
34: .vb
35: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
36: PetscLimiterSetType(PetscLimiter, "my_lim");
37: .ve
38: or at runtime via the option
39: .vb
40: -petsclimiter_type my_lim
41: .ve
43: Level: advanced
45: Note:
46: `PetscLimiterRegister()` may be called multiple times to add several user-defined PetscLimiters
48: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterRegisterAll()`, `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: /*@
58: PetscLimiterSetType - Builds a `PetscLimiter` for a given `PetscLimiterType`
60: Collective
62: Input Parameters:
63: + lim - The `PetscLimiter` object
64: - name - The kind of limiter
66: Options Database Key:
67: . -petsclimiter_type type - Sets the PetscLimiter type; use -help for a list of available types
69: Level: intermediate
71: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
72: @*/
73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
74: {
75: PetscErrorCode (*r)(PetscLimiter);
76: PetscBool match;
78: PetscFunctionBegin;
80: PetscCall(PetscObjectTypeCompare((PetscObject)lim, name, &match));
81: if (match) PetscFunctionReturn(PETSC_SUCCESS);
83: PetscCall(PetscLimiterRegisterAll());
84: PetscCall(PetscFunctionListFind(PetscLimiterList, name, &r));
85: PetscCheck(r, PetscObjectComm((PetscObject)lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
87: PetscTryTypeMethod(lim, destroy);
88: lim->ops->destroy = NULL;
90: PetscCall((*r)(lim));
91: PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@
96: PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.
98: Not Collective
100: Input Parameter:
101: . lim - The `PetscLimiter`
103: Output Parameter:
104: . name - The `PetscLimiterType`
106: Level: intermediate
108: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
112: PetscFunctionBegin;
114: PetscAssertPointer(name, 2);
115: PetscCall(PetscLimiterRegisterAll());
116: *name = ((PetscObject)lim)->type_name;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database
123: Collective
125: Input Parameters:
126: + A - the `PetscLimiter` object to view
127: . obj - Optional object that provides the options prefix to use
128: - name - command line option name
130: Options Database Key:
131: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
133: Level: intermediate
135: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
136: @*/
137: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
138: {
139: PetscFunctionBegin;
141: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@
146: PetscLimiterView - Views a `PetscLimiter`
148: Collective
150: Input Parameters:
151: + lim - the `PetscLimiter` object to view
152: - v - the viewer
154: Level: beginner
156: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
157: @*/
158: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
159: {
160: PetscFunctionBegin;
162: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
163: PetscTryTypeMethod(lim, view, v);
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*@
168: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
170: Collective
172: Input Parameter:
173: . lim - the `PetscLimiter` object to set options for
175: Level: intermediate
177: .seealso: `PetscLimiter`, `PetscLimiterView()`
178: @*/
179: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
180: {
181: const char *defaultType;
182: char name[256];
183: PetscBool flg;
185: PetscFunctionBegin;
187: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
188: else defaultType = ((PetscObject)lim)->type_name;
189: PetscCall(PetscLimiterRegisterAll());
191: PetscObjectOptionsBegin((PetscObject)lim);
192: PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg));
193: if (flg) {
194: PetscCall(PetscLimiterSetType(lim, name));
195: } else if (!((PetscObject)lim)->type_name) {
196: PetscCall(PetscLimiterSetType(lim, defaultType));
197: }
198: PetscTryTypeMethod(lim, setfromoptions);
199: /* process any options handlers added with PetscObjectAddOptionsHandler() */
200: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
201: PetscOptionsEnd();
202: PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*@
207: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
209: Collective
211: Input Parameter:
212: . lim - the `PetscLimiter` object to setup
214: Level: intermediate
216: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscLimiterDestroy()`
217: @*/
218: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
219: {
220: PetscFunctionBegin;
222: PetscTryTypeMethod(lim, setup);
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@
227: PetscLimiterDestroy - Destroys a `PetscLimiter` object
229: Collective
231: Input Parameter:
232: . lim - the `PetscLimiter` object to destroy
234: Level: beginner
236: .seealso: `PetscLimiter`, `PetscLimiterView()`
237: @*/
238: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
239: {
240: PetscFunctionBegin;
241: if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);
244: if (--((PetscObject)*lim)->refct > 0) {
245: *lim = NULL;
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
248: ((PetscObject)*lim)->refct = 0;
250: PetscTryTypeMethod(*lim, destroy);
251: PetscCall(PetscHeaderDestroy(lim));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
258: Collective
260: Input Parameter:
261: . comm - The communicator for the `PetscLimiter` object
263: Output Parameter:
264: . lim - The `PetscLimiter` object
266: Level: beginner
268: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
269: @*/
270: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
271: {
272: PetscLimiter l;
274: PetscFunctionBegin;
275: PetscAssertPointer(lim, 2);
276: PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
277: PetscCall(PetscFVInitializePackage());
279: PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));
281: *lim = l;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: PetscLimiterLimit - Limit the flux
288: Input Parameters:
289: + lim - The `PetscLimiter`
290: - flim - The input field
292: Output Parameter:
293: . phi - The limited field
295: Level: beginner
297: Note:
298: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
299: .vb
300: The classical flux-limited formulation is psi(r) where
302: r = (u[0] - u[-1]) / (u[1] - u[0])
304: The second order TVD region is bounded by
306: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
308: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
309: phi(r)(r+1)/2 in which the reconstructed interface values are
311: u(v) = u[0] + phi(r) (grad u)[0] v
313: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
315: 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))
317: For a nicer symmetric formulation, rewrite in terms of
319: f = (u[0] - u[-1]) / (u[1] - u[-1])
321: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
323: phi(r) = phi(1/r)
325: becomes
327: w(f) = w(1-f).
329: The limiters below implement this final form w(f). The reference methods are
331: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
332: .ve
334: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
335: @*/
336: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
337: {
338: PetscFunctionBegin;
340: PetscAssertPointer(phi, 3);
341: PetscUseTypeMethod(lim, limit, flim, phi);
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
346: {
347: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
349: PetscFunctionBegin;
350: PetscCall(PetscFree(l));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
355: {
356: PetscViewerFormat format;
358: PetscFunctionBegin;
359: PetscCall(PetscViewerGetFormat(viewer, &format));
360: PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
365: {
366: PetscBool isascii;
368: PetscFunctionBegin;
371: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
372: if (isascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
377: {
378: PetscFunctionBegin;
379: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
384: {
385: PetscFunctionBegin;
386: lim->ops->view = PetscLimiterView_Sin;
387: lim->ops->destroy = PetscLimiterDestroy_Sin;
388: lim->ops->limit = PetscLimiterLimit_Sin;
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*MC
393: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
395: Level: intermediate
397: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
398: M*/
400: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
401: {
402: PetscLimiter_Sin *l;
404: PetscFunctionBegin;
406: PetscCall(PetscNew(&l));
407: lim->data = l;
409: PetscCall(PetscLimiterInitialize_Sin(lim));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
414: {
415: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
417: PetscFunctionBegin;
418: PetscCall(PetscFree(l));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
423: {
424: PetscViewerFormat format;
426: PetscFunctionBegin;
427: PetscCall(PetscViewerGetFormat(viewer, &format));
428: PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
433: {
434: PetscBool isascii;
436: PetscFunctionBegin;
439: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
440: if (isascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
445: {
446: PetscFunctionBegin;
447: *phi = 0.0;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
452: {
453: PetscFunctionBegin;
454: lim->ops->view = PetscLimiterView_Zero;
455: lim->ops->destroy = PetscLimiterDestroy_Zero;
456: lim->ops->limit = PetscLimiterLimit_Zero;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*MC
461: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
463: Level: intermediate
465: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
466: M*/
468: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
469: {
470: PetscLimiter_Zero *l;
472: PetscFunctionBegin;
474: PetscCall(PetscNew(&l));
475: lim->data = l;
477: PetscCall(PetscLimiterInitialize_Zero(lim));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
482: {
483: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
485: PetscFunctionBegin;
486: PetscCall(PetscFree(l));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
491: {
492: PetscViewerFormat format;
494: PetscFunctionBegin;
495: PetscCall(PetscViewerGetFormat(viewer, &format));
496: PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
501: {
502: PetscBool isascii;
504: PetscFunctionBegin;
507: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
508: if (isascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
513: {
514: PetscFunctionBegin;
515: *phi = 1.0;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
520: {
521: PetscFunctionBegin;
522: lim->ops->view = PetscLimiterView_None;
523: lim->ops->destroy = PetscLimiterDestroy_None;
524: lim->ops->limit = PetscLimiterLimit_None;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: /*MC
529: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
531: Level: intermediate
533: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
534: M*/
536: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
537: {
538: PetscLimiter_None *l;
540: PetscFunctionBegin;
542: PetscCall(PetscNew(&l));
543: lim->data = l;
545: PetscCall(PetscLimiterInitialize_None(lim));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
550: {
551: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
553: PetscFunctionBegin;
554: PetscCall(PetscFree(l));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
559: {
560: PetscViewerFormat format;
562: PetscFunctionBegin;
563: PetscCall(PetscViewerGetFormat(viewer, &format));
564: PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
569: {
570: PetscBool isascii;
572: PetscFunctionBegin;
575: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
576: if (isascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
581: {
582: PetscFunctionBegin;
583: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
588: {
589: PetscFunctionBegin;
590: lim->ops->view = PetscLimiterView_Minmod;
591: lim->ops->destroy = PetscLimiterDestroy_Minmod;
592: lim->ops->limit = PetscLimiterLimit_Minmod;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*MC
597: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
599: Level: intermediate
601: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
602: M*/
604: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
605: {
606: PetscLimiter_Minmod *l;
608: PetscFunctionBegin;
610: PetscCall(PetscNew(&l));
611: lim->data = l;
613: PetscCall(PetscLimiterInitialize_Minmod(lim));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
618: {
619: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
621: PetscFunctionBegin;
622: PetscCall(PetscFree(l));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
627: {
628: PetscViewerFormat format;
630: PetscFunctionBegin;
631: PetscCall(PetscViewerGetFormat(viewer, &format));
632: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
637: {
638: PetscBool isascii;
640: PetscFunctionBegin;
643: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
644: if (isascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
649: {
650: PetscFunctionBegin;
651: *phi = PetscMax(0, 4 * f * (1 - f));
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
656: {
657: PetscFunctionBegin;
658: lim->ops->view = PetscLimiterView_VanLeer;
659: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
660: lim->ops->limit = PetscLimiterLimit_VanLeer;
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*MC
665: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
667: Level: intermediate
669: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
670: M*/
672: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
673: {
674: PetscLimiter_VanLeer *l;
676: PetscFunctionBegin;
678: PetscCall(PetscNew(&l));
679: lim->data = l;
681: PetscCall(PetscLimiterInitialize_VanLeer(lim));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
686: {
687: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
689: PetscFunctionBegin;
690: PetscCall(PetscFree(l));
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
695: {
696: PetscViewerFormat format;
698: PetscFunctionBegin;
699: PetscCall(PetscViewerGetFormat(viewer, &format));
700: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
705: {
706: PetscBool isascii;
708: PetscFunctionBegin;
711: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
712: if (isascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
717: {
718: PetscFunctionBegin;
719: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
724: {
725: PetscFunctionBegin;
726: lim->ops->view = PetscLimiterView_VanAlbada;
727: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
728: lim->ops->limit = PetscLimiterLimit_VanAlbada;
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: /*MC
733: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
735: Level: intermediate
737: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
738: M*/
740: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
741: {
742: PetscLimiter_VanAlbada *l;
744: PetscFunctionBegin;
746: PetscCall(PetscNew(&l));
747: lim->data = l;
749: PetscCall(PetscLimiterInitialize_VanAlbada(lim));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
754: {
755: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
757: PetscFunctionBegin;
758: PetscCall(PetscFree(l));
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
763: {
764: PetscViewerFormat format;
766: PetscFunctionBegin;
767: PetscCall(PetscViewerGetFormat(viewer, &format));
768: PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
773: {
774: PetscBool isascii;
776: PetscFunctionBegin;
779: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
780: if (isascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
785: {
786: PetscFunctionBegin;
787: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
792: {
793: PetscFunctionBegin;
794: lim->ops->view = PetscLimiterView_Superbee;
795: lim->ops->destroy = PetscLimiterDestroy_Superbee;
796: lim->ops->limit = PetscLimiterLimit_Superbee;
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: /*MC
801: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
803: Level: intermediate
805: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
806: M*/
808: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
809: {
810: PetscLimiter_Superbee *l;
812: PetscFunctionBegin;
814: PetscCall(PetscNew(&l));
815: lim->data = l;
817: PetscCall(PetscLimiterInitialize_Superbee(lim));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
822: {
823: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
825: PetscFunctionBegin;
826: PetscCall(PetscFree(l));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
831: {
832: PetscViewerFormat format;
834: PetscFunctionBegin;
835: PetscCall(PetscViewerGetFormat(viewer, &format));
836: PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
841: {
842: PetscBool isascii;
844: PetscFunctionBegin;
847: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
848: if (isascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /* aka Barth-Jespersen */
853: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
854: {
855: PetscFunctionBegin;
856: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
861: {
862: PetscFunctionBegin;
863: lim->ops->view = PetscLimiterView_MC;
864: lim->ops->destroy = PetscLimiterDestroy_MC;
865: lim->ops->limit = PetscLimiterLimit_MC;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*MC
870: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
872: Level: intermediate
874: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
875: M*/
877: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
878: {
879: PetscLimiter_MC *l;
881: PetscFunctionBegin;
883: PetscCall(PetscNew(&l));
884: lim->data = l;
886: PetscCall(PetscLimiterInitialize_MC(lim));
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: PetscClassId PETSCFV_CLASSID = 0;
892: PetscFunctionList PetscFVList = NULL;
893: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
895: /*@C
896: PetscFVRegister - Adds a new `PetscFV` implementation
898: Not Collective, No Fortran Support
900: Input Parameters:
901: + sname - The name of a new user-defined creation routine
902: - function - The creation routine itself
904: Example Usage:
905: .vb
906: PetscFVRegister("my_fv", MyPetscFVCreate);
907: .ve
909: Then, your PetscFV type can be chosen with the procedural interface via
910: .vb
911: PetscFVCreate(MPI_Comm, PetscFV *);
912: PetscFVSetType(PetscFV, "my_fv");
913: .ve
914: or at runtime via the option
915: .vb
916: -petscfv_type my_fv
917: .ve
919: Level: advanced
921: Note:
922: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
924: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
925: @*/
926: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
927: {
928: PetscFunctionBegin;
929: PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: /*@
934: PetscFVSetType - Builds a particular `PetscFV`
936: Collective
938: Input Parameters:
939: + fvm - The `PetscFV` object
940: - name - The type of FVM space
942: Options Database Key:
943: . -petscfv_type type - Sets the `PetscFVType`; use -help for a list of available types
945: Level: intermediate
947: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
948: @*/
949: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
950: {
951: PetscErrorCode (*r)(PetscFV);
952: PetscBool match;
954: PetscFunctionBegin;
956: PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
957: if (match) PetscFunctionReturn(PETSC_SUCCESS);
959: PetscCall(PetscFVRegisterAll());
960: PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
961: PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
963: PetscTryTypeMethod(fvm, destroy);
964: fvm->ops->destroy = NULL;
966: PetscCall((*r)(fvm));
967: PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
974: Not Collective
976: Input Parameter:
977: . fvm - The `PetscFV`
979: Output Parameter:
980: . name - The `PetscFVType` name
982: Level: intermediate
984: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
985: @*/
986: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
987: {
988: PetscFunctionBegin;
990: PetscAssertPointer(name, 2);
991: PetscCall(PetscFVRegisterAll());
992: *name = ((PetscObject)fvm)->type_name;
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@
997: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
999: Collective
1001: Input Parameters:
1002: + A - the `PetscFV` object
1003: . obj - Optional object that provides the options prefix
1004: - name - command line option name
1006: Options Database Key:
1007: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
1009: Level: intermediate
1011: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
1012: @*/
1013: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1014: {
1015: PetscFunctionBegin;
1017: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: /*@
1022: PetscFVView - Views a `PetscFV`
1024: Collective
1026: Input Parameters:
1027: + fvm - the `PetscFV` object to view
1028: - v - the viewer
1030: Level: beginner
1032: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1033: @*/
1034: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1035: {
1036: PetscFunctionBegin;
1038: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1039: PetscTryTypeMethod(fvm, view, v);
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: /*@
1044: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
1046: Collective
1048: Input Parameter:
1049: . fvm - the `PetscFV` object to set options for
1051: Options Database Key:
1052: . -petscfv_compute_gradients (true|false) - Determines whether cell gradients are calculated
1054: Level: intermediate
1056: .seealso: `PetscFV`, `PetscFVView()`
1057: @*/
1058: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1059: {
1060: const char *defaultType;
1061: char name[256];
1062: PetscBool flg;
1064: PetscFunctionBegin;
1066: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1067: else defaultType = ((PetscObject)fvm)->type_name;
1068: PetscCall(PetscFVRegisterAll());
1070: PetscObjectOptionsBegin((PetscObject)fvm);
1071: PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg));
1072: if (flg) PetscCall(PetscFVSetType(fvm, name));
1073: else if (!((PetscObject)fvm)->type_name) PetscCall(PetscFVSetType(fvm, defaultType));
1074: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1075: PetscTryTypeMethod(fvm, setfromoptions);
1076: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1077: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1078: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1079: PetscOptionsEnd();
1080: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /*@
1085: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1087: Collective
1089: Input Parameter:
1090: . fvm - the `PetscFV` object to setup
1092: Level: intermediate
1094: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1095: @*/
1096: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1097: {
1098: PetscFunctionBegin;
1100: PetscCall(PetscLimiterSetUp(fvm->limiter));
1101: PetscTryTypeMethod(fvm, setup);
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: /*@
1106: PetscFVDestroy - Destroys a `PetscFV` object
1108: Collective
1110: Input Parameter:
1111: . fvm - the `PetscFV` object to destroy
1113: Level: beginner
1115: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1116: @*/
1117: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1118: {
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 (PetscInt 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: PetscCall(PetscFVInitializePackage());
1165: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1166: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1167: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1168: f->numComponents = 1;
1169: f->dim = 0;
1170: f->computeGradients = PETSC_FALSE;
1171: f->fluxWork = NULL;
1172: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1174: *fvm = f;
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: /*@
1179: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1181: Logically Collective
1183: Input Parameters:
1184: + fvm - the `PetscFV` object
1185: - lim - The `PetscLimiter`
1187: Level: intermediate
1189: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1190: @*/
1191: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1192: {
1193: PetscFunctionBegin;
1196: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1197: PetscCall(PetscObjectReference((PetscObject)lim));
1198: fvm->limiter = lim;
1199: PetscFunctionReturn(PETSC_SUCCESS);
1200: }
1202: /*@
1203: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1205: Not Collective
1207: Input Parameter:
1208: . fvm - the `PetscFV` object
1210: Output Parameter:
1211: . lim - The `PetscLimiter`
1213: Level: intermediate
1215: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1216: @*/
1217: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1218: {
1219: PetscFunctionBegin;
1221: PetscAssertPointer(lim, 2);
1222: *lim = fvm->limiter;
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: /*@
1227: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1229: Logically Collective
1231: Input Parameters:
1232: + fvm - the `PetscFV` object
1233: - comp - The number of components
1235: Level: intermediate
1237: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1238: @*/
1239: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1240: {
1241: PetscFunctionBegin;
1243: if (fvm->numComponents != comp) {
1244: for (PetscInt i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1245: PetscCall(PetscFree(fvm->componentNames));
1246: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1247: }
1248: fvm->numComponents = comp;
1249: PetscCall(PetscFree(fvm->fluxWork));
1250: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: /*@
1255: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1257: Not Collective
1259: Input Parameter:
1260: . fvm - the `PetscFV` object
1262: Output Parameter:
1263: . comp - The number of components
1265: Level: intermediate
1267: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1268: @*/
1269: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1270: {
1271: PetscFunctionBegin;
1273: PetscAssertPointer(comp, 2);
1274: *comp = fvm->numComponents;
1275: PetscFunctionReturn(PETSC_SUCCESS);
1276: }
1278: /*@
1279: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1281: Logically Collective
1283: Input Parameters:
1284: + fvm - the `PetscFV` object
1285: . comp - the component number
1286: - name - the component name
1288: Level: intermediate
1290: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1291: @*/
1292: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1293: {
1294: PetscFunctionBegin;
1295: PetscCall(PetscFree(fvm->componentNames[comp]));
1296: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1297: PetscFunctionReturn(PETSC_SUCCESS);
1298: }
1300: /*@
1301: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1303: Logically Collective
1305: Input Parameters:
1306: + fvm - the `PetscFV` object
1307: - comp - the component number
1309: Output Parameter:
1310: . name - the component name
1312: Level: intermediate
1314: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1315: @*/
1316: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char *name[])
1317: {
1318: PetscFunctionBegin;
1319: *name = fvm->componentNames[comp];
1320: PetscFunctionReturn(PETSC_SUCCESS);
1321: }
1323: /*@
1324: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1326: Logically Collective
1328: Input Parameters:
1329: + fvm - the `PetscFV` object
1330: - dim - The spatial dimension
1332: Level: intermediate
1334: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
1335: @*/
1336: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1337: {
1338: PetscFunctionBegin;
1340: fvm->dim = dim;
1341: PetscFunctionReturn(PETSC_SUCCESS);
1342: }
1344: /*@
1345: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1347: Not Collective
1349: Input Parameter:
1350: . fvm - the `PetscFV` object
1352: Output Parameter:
1353: . dim - The spatial dimension
1355: Level: intermediate
1357: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1358: @*/
1359: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1360: {
1361: PetscFunctionBegin;
1363: PetscAssertPointer(dim, 2);
1364: *dim = fvm->dim;
1365: PetscFunctionReturn(PETSC_SUCCESS);
1366: }
1368: /*@
1369: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1371: Logically Collective
1373: Input Parameters:
1374: + fvm - the `PetscFV` object
1375: - computeGradients - Flag to compute cell gradients
1377: Level: intermediate
1379: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1380: @*/
1381: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1382: {
1383: PetscFunctionBegin;
1385: fvm->computeGradients = computeGradients;
1386: PetscFunctionReturn(PETSC_SUCCESS);
1387: }
1389: /*@
1390: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1392: Not Collective
1394: Input Parameter:
1395: . fvm - the `PetscFV` object
1397: Output Parameter:
1398: . computeGradients - Flag to compute cell gradients
1400: Level: intermediate
1402: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1403: @*/
1404: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1405: {
1406: PetscFunctionBegin;
1408: PetscAssertPointer(computeGradients, 2);
1409: *computeGradients = fvm->computeGradients;
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: /*@
1414: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1416: Logically Collective
1418: Input Parameters:
1419: + fvm - the `PetscFV` object
1420: - q - The `PetscQuadrature`
1422: Level: intermediate
1424: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1425: @*/
1426: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1427: {
1428: PetscFunctionBegin;
1430: PetscCall(PetscObjectReference((PetscObject)q));
1431: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1432: fvm->quadrature = q;
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1439: Not Collective
1441: Input Parameter:
1442: . fvm - the `PetscFV` object
1444: Output Parameter:
1445: . q - The `PetscQuadrature`
1447: Level: intermediate
1449: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1450: @*/
1451: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1452: {
1453: PetscFunctionBegin;
1455: PetscAssertPointer(q, 2);
1456: if (!fvm->quadrature) {
1457: /* Create default 1-point quadrature */
1458: PetscReal *points, *weights;
1460: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1461: PetscCall(PetscCalloc1(fvm->dim, &points));
1462: PetscCall(PetscMalloc1(1, &weights));
1463: weights[0] = 1.0;
1464: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1465: }
1466: *q = fvm->quadrature;
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }
1470: /*@
1471: PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`
1473: Not Collective
1475: Input Parameters:
1476: + fvm - The `PetscFV` object
1477: - ct - The `DMPolytopeType` for the cell
1479: Level: intermediate
1481: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1482: @*/
1483: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1484: {
1485: DM K;
1486: PetscInt dim, Nc;
1488: PetscFunctionBegin;
1489: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1490: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1491: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1492: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1493: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1494: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1495: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1496: PetscCall(DMDestroy(&K));
1497: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1498: // Should we be using PetscFVGetQuadrature() here?
1499: for (PetscInt c = 0; c < Nc; ++c) {
1500: PetscQuadrature qc;
1501: PetscReal *points, *weights;
1503: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1504: PetscCall(PetscCalloc1(dim, &points));
1505: PetscCall(PetscCalloc1(Nc, &weights));
1506: weights[c] = 1.0;
1507: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1508: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1509: PetscCall(PetscQuadratureDestroy(&qc));
1510: }
1511: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1512: PetscFunctionReturn(PETSC_SUCCESS);
1513: }
1515: /*@
1516: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1518: Not Collective
1520: Input Parameter:
1521: . fvm - The `PetscFV` object
1523: Output Parameter:
1524: . sp - The `PetscDualSpace` object
1526: Level: intermediate
1528: Developer Notes:
1529: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1531: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1532: @*/
1533: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1534: {
1535: PetscFunctionBegin;
1537: PetscAssertPointer(sp, 2);
1538: if (!fvm->dualSpace) {
1539: PetscInt dim;
1541: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1542: PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1543: }
1544: *sp = fvm->dualSpace;
1545: PetscFunctionReturn(PETSC_SUCCESS);
1546: }
1548: /*@
1549: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1551: Not Collective
1553: Input Parameters:
1554: + fvm - The `PetscFV` object
1555: - sp - The `PetscDualSpace` object
1557: Level: intermediate
1559: Note:
1560: A simple dual space is provided automatically, and the user typically will not need to override it.
1562: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1563: @*/
1564: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1565: {
1566: PetscFunctionBegin;
1569: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1570: fvm->dualSpace = sp;
1571: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1572: PetscFunctionReturn(PETSC_SUCCESS);
1573: }
1575: /*@C
1576: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1578: Not Collective
1580: Input Parameter:
1581: . fvm - The `PetscFV` object
1583: Output Parameter:
1584: . T - The basis function values and derivatives at quadrature points
1586: Level: intermediate
1588: Note:
1589: .vb
1590: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1591: 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
1592: 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
1593: .ve
1595: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1596: @*/
1597: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1598: {
1599: PetscInt npoints;
1600: const PetscReal *points;
1602: PetscFunctionBegin;
1604: PetscAssertPointer(T, 2);
1605: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1606: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1607: *T = fvm->T;
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: /*@C
1612: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1614: Not Collective
1616: Input Parameters:
1617: + fvm - The `PetscFV` object
1618: . nrepl - The number of replicas
1619: . npoints - The number of tabulation points in a replica
1620: . points - The tabulation point coordinates
1621: - K - The order of derivative to tabulate
1623: Output Parameter:
1624: . T - The basis function values and derivative at tabulation points
1626: Level: intermediate
1628: Note:
1629: .vb
1630: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1631: 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
1632: 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
1633: .ve
1635: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1636: @*/
1637: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1638: {
1639: PetscInt pdim; // Dimension of approximation space P
1640: PetscInt cdim; // Spatial dimension
1641: PetscInt Nc; // Field components
1642: PetscInt k, p, d, c, e;
1644: PetscFunctionBegin;
1645: if (!npoints || K < 0) {
1646: *T = NULL;
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: PetscAssertPointer(points, 4);
1651: PetscAssertPointer(T, 6);
1652: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1653: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1654: pdim = Nc;
1655: PetscCall(PetscMalloc1(1, T));
1656: (*T)->K = !cdim ? 0 : K;
1657: (*T)->Nr = nrepl;
1658: (*T)->Np = npoints;
1659: (*T)->Nb = pdim;
1660: (*T)->Nc = Nc;
1661: (*T)->cdim = cdim;
1662: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1663: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1664: if (K >= 0) {
1665: for (p = 0; p < nrepl * npoints; ++p)
1666: for (d = 0; d < pdim; ++d)
1667: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1668: }
1669: if (K >= 1) {
1670: for (p = 0; p < nrepl * npoints; ++p)
1671: for (d = 0; d < pdim; ++d)
1672: for (c = 0; c < Nc; ++c)
1673: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1674: }
1675: if (K >= 2) {
1676: for (p = 0; p < nrepl * npoints; ++p)
1677: for (d = 0; d < pdim; ++d)
1678: for (c = 0; c < Nc; ++c)
1679: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1680: }
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: /*@
1685: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1687: Input Parameters:
1688: + fvm - The `PetscFV` object
1689: . numFaces - The number of cell faces which are not constrained
1690: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1692: Output Parameter:
1693: . grad - the gradient
1695: Level: advanced
1697: .seealso: `PetscFV`, `PetscFVCreate()`
1698: @*/
1699: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1700: {
1701: PetscFunctionBegin;
1703: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /*@C
1708: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1710: Not Collective
1712: Input Parameters:
1713: + fvm - The `PetscFV` object for the field being integrated
1714: . prob - The `PetscDS` specifying the discretizations and continuum functions
1715: . field - The field being integrated
1716: . Nf - The number of faces in the chunk
1717: . fgeom - The face geometry for each face in the chunk
1718: . neighborVol - The volume for each pair of cells in the chunk
1719: . uL - The state from the cell on the left
1720: - uR - The state from the cell on the right
1722: Output Parameters:
1723: + fluxL - the left fluxes for each face
1724: - fluxR - the right fluxes for each face
1726: Level: developer
1728: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1729: @*/
1730: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1731: {
1732: PetscFunctionBegin;
1734: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1735: PetscFunctionReturn(PETSC_SUCCESS);
1736: }
1738: /*@
1739: PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.
1741: Input Parameter:
1742: . fv - The initial `PetscFV`
1744: Output Parameter:
1745: . fvNew - A clone of the `PetscFV`
1747: Level: advanced
1749: Notes:
1750: This is typically used to change the number of components.
1752: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1753: @*/
1754: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1755: {
1756: PetscDualSpace Q;
1757: DM K;
1758: PetscQuadrature q;
1759: PetscInt Nc, cdim;
1761: PetscFunctionBegin;
1762: PetscCall(PetscFVGetDualSpace(fv, &Q));
1763: PetscCall(PetscFVGetQuadrature(fv, &q));
1764: PetscCall(PetscDualSpaceGetDM(Q, &K));
1766: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1767: PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1768: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1769: PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1770: PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1771: PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1772: PetscCall(PetscFVSetQuadrature(*fvNew, q));
1773: PetscFunctionReturn(PETSC_SUCCESS);
1774: }
1776: /*@
1777: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1778: smaller copies.
1780: Input Parameter:
1781: . fv - The initial `PetscFV`
1783: Output Parameter:
1784: . fvRef - The refined `PetscFV`
1786: Level: advanced
1788: Notes:
1789: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1790: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1791: interpolation between regularly refined meshes.
1793: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1794: @*/
1795: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1796: {
1797: PetscDualSpace Q, Qref;
1798: DM K, Kref;
1799: PetscQuadrature q, qref;
1800: DMPolytopeType ct;
1801: DMPlexTransform tr;
1802: PetscReal *v0;
1803: PetscReal *jac, *invjac;
1804: PetscInt numComp, numSubelements, s;
1806: PetscFunctionBegin;
1807: PetscCall(PetscFVGetDualSpace(fv, &Q));
1808: PetscCall(PetscFVGetQuadrature(fv, &q));
1809: PetscCall(PetscDualSpaceGetDM(Q, &K));
1810: /* Create dual space */
1811: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1812: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1813: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1814: PetscCall(DMDestroy(&Kref));
1815: PetscCall(PetscDualSpaceSetUp(Qref));
1816: /* Create volume */
1817: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1818: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1819: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1820: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1821: PetscCall(PetscFVSetUp(*fvRef));
1822: /* Create quadrature */
1823: PetscCall(DMPlexGetCellType(K, 0, &ct));
1824: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1825: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1826: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1827: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1828: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1829: for (s = 0; s < numSubelements; ++s) {
1830: PetscQuadrature qs;
1831: const PetscReal *points, *weights;
1832: PetscReal *p, *w;
1833: PetscInt dim, Nc, npoints, np;
1835: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1836: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1837: np = npoints / numSubelements;
1838: PetscCall(PetscMalloc1(np * dim, &p));
1839: PetscCall(PetscMalloc1(np * Nc, &w));
1840: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1841: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1842: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1843: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1844: PetscCall(PetscQuadratureDestroy(&qs));
1845: }
1846: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1847: PetscCall(DMPlexTransformDestroy(&tr));
1848: PetscCall(PetscQuadratureDestroy(&qref));
1849: PetscCall(PetscDualSpaceDestroy(&Qref));
1850: PetscFunctionReturn(PETSC_SUCCESS);
1851: }
1853: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1854: {
1855: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1857: PetscFunctionBegin;
1858: PetscCall(PetscFree(b));
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1863: {
1864: PetscInt Nc;
1865: PetscViewerFormat format;
1867: PetscFunctionBegin;
1868: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1869: PetscCall(PetscViewerGetFormat(viewer, &format));
1870: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1871: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1872: for (PetscInt c = 0; c < Nc; c++) {
1873: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1874: }
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1879: {
1880: PetscBool isascii;
1882: PetscFunctionBegin;
1885: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1886: if (isascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1887: PetscFunctionReturn(PETSC_SUCCESS);
1888: }
1890: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1891: {
1892: PetscInt dim;
1894: PetscFunctionBegin;
1895: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1896: for (PetscInt f = 0; f < numFaces; ++f) {
1897: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1898: }
1899: PetscFunctionReturn(PETSC_SUCCESS);
1900: }
1902: /*
1903: neighborVol[f*2+0] contains the left geom
1904: neighborVol[f*2+1] contains the right geom
1905: */
1906: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1907: {
1908: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1909: void *rctx;
1910: PetscScalar *flux = fvm->fluxWork;
1911: const PetscScalar *constants;
1912: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1914: PetscFunctionBegin;
1915: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1916: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1917: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1918: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1919: PetscCall(PetscDSGetContext(prob, field, &rctx));
1920: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1921: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1922: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1923: for (f = 0; f < Nf; ++f) {
1924: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1925: for (d = 0; d < pdim; ++d) {
1926: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1927: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1928: }
1929: }
1930: PetscFunctionReturn(PETSC_SUCCESS);
1931: }
1933: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1934: {
1935: PetscFunctionBegin;
1936: fvm->ops->setfromoptions = NULL;
1937: fvm->ops->view = PetscFVView_Upwind;
1938: fvm->ops->destroy = PetscFVDestroy_Upwind;
1939: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1940: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1941: PetscFunctionReturn(PETSC_SUCCESS);
1942: }
1944: /*MC
1945: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1947: Level: intermediate
1949: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1950: M*/
1952: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1953: {
1954: PetscFV_Upwind *b;
1956: PetscFunctionBegin;
1958: PetscCall(PetscNew(&b));
1959: fvm->data = b;
1961: PetscCall(PetscFVInitialize_Upwind(fvm));
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: #include <petscblaslapack.h>
1967: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1968: {
1969: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1971: PetscFunctionBegin;
1972: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1973: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1974: PetscCall(PetscFree(ls));
1975: PetscFunctionReturn(PETSC_SUCCESS);
1976: }
1978: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1979: {
1980: PetscInt Nc;
1981: PetscViewerFormat format;
1983: PetscFunctionBegin;
1984: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1985: PetscCall(PetscViewerGetFormat(viewer, &format));
1986: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1987: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1988: for (PetscInt c = 0; c < Nc; c++) {
1989: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1990: }
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1995: {
1996: PetscBool isascii;
1998: PetscFunctionBegin;
2001: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2002: if (isascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2003: PetscFunctionReturn(PETSC_SUCCESS);
2004: }
2006: /* Overwrites A. Can only handle full-rank problems with m>=n */
2007: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2008: {
2009: PetscBool debug = PETSC_FALSE;
2010: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2011: PetscScalar *R, *Q, *Aback, Alpha;
2013: PetscFunctionBegin;
2014: if (debug) {
2015: PetscCall(PetscMalloc1(m * n, &Aback));
2016: PetscCall(PetscArraycpy(Aback, A, m * n));
2017: }
2019: PetscCall(PetscBLASIntCast(m, &M));
2020: PetscCall(PetscBLASIntCast(n, &N));
2021: PetscCall(PetscBLASIntCast(mstride, &lda));
2022: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2023: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2024: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2025: PetscCall(PetscFPTrapPop());
2026: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2027: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2029: /* Extract an explicit representation of Q */
2030: Q = Ainv;
2031: PetscCall(PetscArraycpy(Q, A, mstride * n));
2032: K = N; /* full rank */
2033: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2034: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2036: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2037: Alpha = 1.0;
2038: ldb = lda;
2039: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2040: /* Ainv is Q, overwritten with inverse */
2042: if (debug) { /* Check that pseudo-inverse worked */
2043: PetscScalar Beta = 0.0;
2044: PetscBLASInt ldc;
2045: K = N;
2046: ldc = N;
2047: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2048: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2049: PetscCall(PetscFree(Aback));
2050: }
2051: PetscFunctionReturn(PETSC_SUCCESS);
2052: }
2054: /* Overwrites A. Can handle degenerate problems and m<n. */
2055: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2056: {
2057: PetscScalar *Brhs;
2058: PetscScalar *tmpwork;
2059: PetscReal rcond;
2060: #if defined(PETSC_USE_COMPLEX)
2061: PetscInt rworkSize;
2062: PetscReal *rwork, *rtau;
2063: #endif
2064: PetscInt i, j, maxmn;
2065: PetscBLASInt M, N, lda, ldb, ldwork;
2066: PetscBLASInt nrhs, irank, info;
2068: PetscFunctionBegin;
2069: /* initialize to identity */
2070: tmpwork = work;
2071: Brhs = Ainv;
2072: maxmn = PetscMax(m, n);
2073: for (j = 0; j < maxmn; j++) {
2074: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2075: }
2077: PetscCall(PetscBLASIntCast(m, &M));
2078: PetscCall(PetscBLASIntCast(n, &N));
2079: PetscCall(PetscBLASIntCast(mstride, &lda));
2080: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2081: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2082: rcond = -1;
2083: nrhs = M;
2084: #if defined(PETSC_USE_COMPLEX)
2085: rworkSize = 5 * PetscMin(M, N);
2086: PetscCall(PetscMalloc1(rworkSize, &rwork));
2087: PetscCall(PetscMalloc1(PetscMin(M, N), &rtau));
2088: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2089: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, rtau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2090: PetscCall(PetscFPTrapPop());
2091: PetscCall(PetscFree(rwork));
2092: for (i = 0; i < PetscMin(M, N); i++) tau[i] = rtau[i];
2093: PetscCall(PetscFree(rtau));
2094: #else
2095: nrhs = M;
2096: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2097: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, tau, &rcond, &irank, tmpwork, &ldwork, &info));
2098: PetscCall(PetscFPTrapPop());
2099: #endif
2100: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2101: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2102: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two collinear points");
2103: PetscFunctionReturn(PETSC_SUCCESS);
2104: }
2106: #if 0
2107: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2108: {
2109: PetscReal grad[2] = {0, 0};
2110: const PetscInt *faces;
2111: PetscInt numFaces, f;
2113: PetscFunctionBegin;
2114: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2115: PetscCall(DMPlexGetCone(dm, cell, &faces));
2116: for (f = 0; f < numFaces; ++f) {
2117: const PetscInt *fcells;
2118: const CellGeom *cg1;
2119: const FaceGeom *fg;
2121: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2122: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2123: for (i = 0; i < 2; ++i) {
2124: PetscScalar du;
2126: if (fcells[i] == c) continue;
2127: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2128: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2129: grad[0] += fg->grad[!i][0] * du;
2130: grad[1] += fg->grad[!i][1] * du;
2131: }
2132: }
2133: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2134: PetscFunctionReturn(PETSC_SUCCESS);
2135: }
2136: #endif
2138: /*
2139: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2141: Input Parameters:
2142: + fvm - The `PetscFV` object
2143: . numFaces - The number of cell faces which are not constrained
2144: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2146: Level: developer
2148: .seealso: `PetscFV`, `PetscFVCreate()`
2149: */
2150: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2151: {
2152: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2153: const PetscBool useSVD = PETSC_TRUE;
2154: const PetscInt maxFaces = ls->maxFaces;
2155: PetscInt dim, f, d;
2157: PetscFunctionBegin;
2158: if (numFaces > maxFaces) {
2159: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2160: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2161: }
2162: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2163: for (f = 0; f < numFaces; ++f) {
2164: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2165: }
2166: /* Overwrites B with garbage, returns Binv in row-major format */
2167: if (useSVD) {
2168: PetscInt maxmn = PetscMax(numFaces, dim);
2169: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2170: /* Binv shaped in column-major, coldim=maxmn.*/
2171: for (f = 0; f < numFaces; ++f) {
2172: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2173: }
2174: } else {
2175: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2176: /* Binv shaped in row-major, rowdim=maxFaces.*/
2177: for (f = 0; f < numFaces; ++f) {
2178: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2179: }
2180: }
2181: PetscFunctionReturn(PETSC_SUCCESS);
2182: }
2184: /*
2185: neighborVol[f*2+0] contains the left geom
2186: neighborVol[f*2+1] contains the right geom
2187: */
2188: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2189: {
2190: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2191: void *rctx;
2192: PetscScalar *flux = fvm->fluxWork;
2193: const PetscScalar *constants;
2194: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2196: PetscFunctionBegin;
2197: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2198: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2199: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2200: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2201: PetscCall(PetscDSGetContext(prob, field, &rctx));
2202: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2203: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2204: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2205: for (f = 0; f < Nf; ++f) {
2206: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2207: for (d = 0; d < pdim; ++d) {
2208: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2209: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2210: }
2211: }
2212: PetscFunctionReturn(PETSC_SUCCESS);
2213: }
2215: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2216: {
2217: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2218: PetscInt dim, m, n, nrhs, minmn, maxmn;
2220: PetscFunctionBegin;
2222: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2223: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2224: ls->maxFaces = maxFaces;
2225: m = ls->maxFaces;
2226: n = dim;
2227: nrhs = ls->maxFaces;
2228: minmn = PetscMin(m, n);
2229: maxmn = PetscMax(m, n);
2230: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2231: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2232: PetscFunctionReturn(PETSC_SUCCESS);
2233: }
2235: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2236: {
2237: PetscFunctionBegin;
2238: fvm->ops->setfromoptions = NULL;
2239: fvm->ops->view = PetscFVView_LeastSquares;
2240: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2241: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2242: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: /*MC
2247: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2249: Level: intermediate
2251: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2252: M*/
2254: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2255: {
2256: PetscFV_LeastSquares *ls;
2258: PetscFunctionBegin;
2260: PetscCall(PetscNew(&ls));
2261: fvm->data = ls;
2263: ls->maxFaces = -1;
2264: ls->workSize = -1;
2265: ls->B = NULL;
2266: ls->Binv = NULL;
2267: ls->tau = NULL;
2268: ls->work = NULL;
2270: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2271: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2272: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2273: PetscFunctionReturn(PETSC_SUCCESS);
2274: }
2276: /*@
2277: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2279: Not Collective
2281: Input Parameters:
2282: + fvm - The `PetscFV` object
2283: - maxFaces - The maximum number of cell faces
2285: Level: intermediate
2287: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2288: @*/
2289: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2290: {
2291: PetscFunctionBegin;
2293: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2294: PetscFunctionReturn(PETSC_SUCCESS);
2295: }