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: PetscInt i;
1121: PetscFunctionBegin;
1122: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1125: if (--((PetscObject)*fvm)->refct > 0) {
1126: *fvm = NULL;
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1129: ((PetscObject)*fvm)->refct = 0;
1131: for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1132: PetscCall(PetscFree((*fvm)->componentNames));
1133: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1134: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1135: PetscCall(PetscFree((*fvm)->fluxWork));
1136: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1137: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1139: PetscTryTypeMethod(*fvm, destroy);
1140: PetscCall(PetscHeaderDestroy(fvm));
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: /*@
1145: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1147: Collective
1149: Input Parameter:
1150: . comm - The communicator for the `PetscFV` object
1152: Output Parameter:
1153: . fvm - The `PetscFV` object
1155: Level: beginner
1157: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1158: @*/
1159: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1160: {
1161: PetscFV f;
1163: PetscFunctionBegin;
1164: PetscAssertPointer(fvm, 2);
1165: PetscCall(PetscFVInitializePackage());
1167: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1168: 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: /*@
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: /*@
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: /*@
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 isascii;
1886: PetscFunctionBegin;
1889: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1890: if (isascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1891: PetscFunctionReturn(PETSC_SUCCESS);
1892: }
1894: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1895: {
1896: PetscInt dim;
1898: PetscFunctionBegin;
1899: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1900: for (PetscInt f = 0; f < numFaces; ++f) {
1901: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1902: }
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: /*
1907: neighborVol[f*2+0] contains the left geom
1908: neighborVol[f*2+1] contains the right geom
1909: */
1910: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1911: {
1912: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1913: void *rctx;
1914: PetscScalar *flux = fvm->fluxWork;
1915: const PetscScalar *constants;
1916: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1918: PetscFunctionBegin;
1919: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1920: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1921: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1922: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1923: PetscCall(PetscDSGetContext(prob, field, &rctx));
1924: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1925: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1926: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1927: for (f = 0; f < Nf; ++f) {
1928: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1929: for (d = 0; d < pdim; ++d) {
1930: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1931: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1932: }
1933: }
1934: PetscFunctionReturn(PETSC_SUCCESS);
1935: }
1937: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1938: {
1939: PetscFunctionBegin;
1940: fvm->ops->setfromoptions = NULL;
1941: fvm->ops->view = PetscFVView_Upwind;
1942: fvm->ops->destroy = PetscFVDestroy_Upwind;
1943: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1944: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1945: PetscFunctionReturn(PETSC_SUCCESS);
1946: }
1948: /*MC
1949: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1951: Level: intermediate
1953: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1954: M*/
1956: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1957: {
1958: PetscFV_Upwind *b;
1960: PetscFunctionBegin;
1962: PetscCall(PetscNew(&b));
1963: fvm->data = b;
1965: PetscCall(PetscFVInitialize_Upwind(fvm));
1966: PetscFunctionReturn(PETSC_SUCCESS);
1967: }
1969: #include <petscblaslapack.h>
1971: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1972: {
1973: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1975: PetscFunctionBegin;
1976: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1977: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1978: PetscCall(PetscFree(ls));
1979: PetscFunctionReturn(PETSC_SUCCESS);
1980: }
1982: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1983: {
1984: PetscInt Nc, c;
1985: PetscViewerFormat format;
1987: PetscFunctionBegin;
1988: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1989: PetscCall(PetscViewerGetFormat(viewer, &format));
1990: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1991: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1992: for (c = 0; c < Nc; c++) {
1993: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1994: }
1995: PetscFunctionReturn(PETSC_SUCCESS);
1996: }
1998: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1999: {
2000: PetscBool isascii;
2002: PetscFunctionBegin;
2005: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2006: if (isascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2007: PetscFunctionReturn(PETSC_SUCCESS);
2008: }
2010: /* Overwrites A. Can only handle full-rank problems with m>=n */
2011: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2012: {
2013: PetscBool debug = PETSC_FALSE;
2014: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2015: PetscScalar *R, *Q, *Aback, Alpha;
2017: PetscFunctionBegin;
2018: if (debug) {
2019: PetscCall(PetscMalloc1(m * n, &Aback));
2020: PetscCall(PetscArraycpy(Aback, A, m * n));
2021: }
2023: PetscCall(PetscBLASIntCast(m, &M));
2024: PetscCall(PetscBLASIntCast(n, &N));
2025: PetscCall(PetscBLASIntCast(mstride, &lda));
2026: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2027: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2028: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2029: PetscCall(PetscFPTrapPop());
2030: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2031: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2033: /* Extract an explicit representation of Q */
2034: Q = Ainv;
2035: PetscCall(PetscArraycpy(Q, A, mstride * n));
2036: K = N; /* full rank */
2037: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2038: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2040: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2041: Alpha = 1.0;
2042: ldb = lda;
2043: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2044: /* Ainv is Q, overwritten with inverse */
2046: if (debug) { /* Check that pseudo-inverse worked */
2047: PetscScalar Beta = 0.0;
2048: PetscBLASInt ldc;
2049: K = N;
2050: ldc = N;
2051: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2052: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2053: PetscCall(PetscFree(Aback));
2054: }
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: /* Overwrites A. Can handle degenerate problems and m<n. */
2059: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2060: {
2061: PetscScalar *Brhs;
2062: PetscScalar *tmpwork;
2063: PetscReal rcond;
2064: #if defined(PETSC_USE_COMPLEX)
2065: PetscInt rworkSize;
2066: PetscReal *rwork, *rtau;
2067: #endif
2068: PetscInt i, j, maxmn;
2069: PetscBLASInt M, N, lda, ldb, ldwork;
2070: PetscBLASInt nrhs, irank, info;
2072: PetscFunctionBegin;
2073: /* initialize to identity */
2074: tmpwork = work;
2075: Brhs = Ainv;
2076: maxmn = PetscMax(m, n);
2077: for (j = 0; j < maxmn; j++) {
2078: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2079: }
2081: PetscCall(PetscBLASIntCast(m, &M));
2082: PetscCall(PetscBLASIntCast(n, &N));
2083: PetscCall(PetscBLASIntCast(mstride, &lda));
2084: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2085: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2086: rcond = -1;
2087: nrhs = M;
2088: #if defined(PETSC_USE_COMPLEX)
2089: rworkSize = 5 * PetscMin(M, N);
2090: PetscCall(PetscMalloc1(rworkSize, &rwork));
2091: PetscCall(PetscMalloc1(PetscMin(M, N), &rtau));
2092: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2093: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, rtau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2094: PetscCall(PetscFPTrapPop());
2095: PetscCall(PetscFree(rwork));
2096: for (i = 0; i < PetscMin(M, N); i++) tau[i] = rtau[i];
2097: PetscCall(PetscFree(rtau));
2098: #else
2099: nrhs = M;
2100: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2101: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, tau, &rcond, &irank, tmpwork, &ldwork, &info));
2102: PetscCall(PetscFPTrapPop());
2103: #endif
2104: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2105: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2106: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two collinear points");
2107: PetscFunctionReturn(PETSC_SUCCESS);
2108: }
2110: #if 0
2111: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2112: {
2113: PetscReal grad[2] = {0, 0};
2114: const PetscInt *faces;
2115: PetscInt numFaces, f;
2117: PetscFunctionBegin;
2118: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2119: PetscCall(DMPlexGetCone(dm, cell, &faces));
2120: for (f = 0; f < numFaces; ++f) {
2121: const PetscInt *fcells;
2122: const CellGeom *cg1;
2123: const FaceGeom *fg;
2125: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2126: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2127: for (i = 0; i < 2; ++i) {
2128: PetscScalar du;
2130: if (fcells[i] == c) continue;
2131: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2132: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2133: grad[0] += fg->grad[!i][0] * du;
2134: grad[1] += fg->grad[!i][1] * du;
2135: }
2136: }
2137: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2138: PetscFunctionReturn(PETSC_SUCCESS);
2139: }
2140: #endif
2142: /*
2143: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2145: Input Parameters:
2146: + fvm - The `PetscFV` object
2147: . numFaces - The number of cell faces which are not constrained
2148: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2150: Level: developer
2152: .seealso: `PetscFV`, `PetscFVCreate()`
2153: */
2154: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2155: {
2156: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2157: const PetscBool useSVD = PETSC_TRUE;
2158: const PetscInt maxFaces = ls->maxFaces;
2159: PetscInt dim, f, d;
2161: PetscFunctionBegin;
2162: if (numFaces > maxFaces) {
2163: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2164: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2165: }
2166: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2167: for (f = 0; f < numFaces; ++f) {
2168: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2169: }
2170: /* Overwrites B with garbage, returns Binv in row-major format */
2171: if (useSVD) {
2172: PetscInt maxmn = PetscMax(numFaces, dim);
2173: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2174: /* Binv shaped in column-major, coldim=maxmn.*/
2175: for (f = 0; f < numFaces; ++f) {
2176: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2177: }
2178: } else {
2179: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2180: /* Binv shaped in row-major, rowdim=maxFaces.*/
2181: for (f = 0; f < numFaces; ++f) {
2182: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2183: }
2184: }
2185: PetscFunctionReturn(PETSC_SUCCESS);
2186: }
2188: /*
2189: neighborVol[f*2+0] contains the left geom
2190: neighborVol[f*2+1] contains the right geom
2191: */
2192: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2193: {
2194: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2195: void *rctx;
2196: PetscScalar *flux = fvm->fluxWork;
2197: const PetscScalar *constants;
2198: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2200: PetscFunctionBegin;
2201: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2202: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2203: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2204: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2205: PetscCall(PetscDSGetContext(prob, field, &rctx));
2206: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2207: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2208: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2209: for (f = 0; f < Nf; ++f) {
2210: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2211: for (d = 0; d < pdim; ++d) {
2212: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2213: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2214: }
2215: }
2216: PetscFunctionReturn(PETSC_SUCCESS);
2217: }
2219: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2220: {
2221: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2222: PetscInt dim, m, n, nrhs, minmn, maxmn;
2224: PetscFunctionBegin;
2226: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2227: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2228: ls->maxFaces = maxFaces;
2229: m = ls->maxFaces;
2230: n = dim;
2231: nrhs = ls->maxFaces;
2232: minmn = PetscMin(m, n);
2233: maxmn = PetscMax(m, n);
2234: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2235: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2236: PetscFunctionReturn(PETSC_SUCCESS);
2237: }
2239: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2240: {
2241: PetscFunctionBegin;
2242: fvm->ops->setfromoptions = NULL;
2243: fvm->ops->view = PetscFVView_LeastSquares;
2244: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2245: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2246: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2247: PetscFunctionReturn(PETSC_SUCCESS);
2248: }
2250: /*MC
2251: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2253: Level: intermediate
2255: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2256: M*/
2258: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2259: {
2260: PetscFV_LeastSquares *ls;
2262: PetscFunctionBegin;
2264: PetscCall(PetscNew(&ls));
2265: fvm->data = ls;
2267: ls->maxFaces = -1;
2268: ls->workSize = -1;
2269: ls->B = NULL;
2270: ls->Binv = NULL;
2271: ls->tau = NULL;
2272: ls->work = NULL;
2274: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2275: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2276: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2277: PetscFunctionReturn(PETSC_SUCCESS);
2278: }
2280: /*@
2281: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2283: Not Collective
2285: Input Parameters:
2286: + fvm - The `PetscFV` object
2287: - maxFaces - The maximum number of cell faces
2289: Level: intermediate
2291: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2292: @*/
2293: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2294: {
2295: PetscFunctionBegin;
2297: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2298: PetscFunctionReturn(PETSC_SUCCESS);
2299: }