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: Level: intermediate
132: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
133: @*/
134: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
135: {
136: PetscFunctionBegin;
138: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*@
143: PetscLimiterView - Views a `PetscLimiter`
145: Collective
147: Input Parameters:
148: + lim - the `PetscLimiter` object to view
149: - v - the viewer
151: Level: beginner
153: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
154: @*/
155: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
156: {
157: PetscFunctionBegin;
159: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
160: PetscTryTypeMethod(lim, view, v);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database
167: Collective
169: Input Parameter:
170: . lim - the `PetscLimiter` object to set options for
172: Level: intermediate
174: .seealso: `PetscLimiter`, `PetscLimiterView()`
175: @*/
176: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
177: {
178: const char *defaultType;
179: char name[256];
180: PetscBool flg;
182: PetscFunctionBegin;
184: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
185: else defaultType = ((PetscObject)lim)->type_name;
186: PetscCall(PetscLimiterRegisterAll());
188: PetscObjectOptionsBegin((PetscObject)lim);
189: PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg));
190: if (flg) {
191: PetscCall(PetscLimiterSetType(lim, name));
192: } else if (!((PetscObject)lim)->type_name) {
193: PetscCall(PetscLimiterSetType(lim, defaultType));
194: }
195: PetscTryTypeMethod(lim, setfromoptions);
196: /* process any options handlers added with PetscObjectAddOptionsHandler() */
197: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
198: PetscOptionsEnd();
199: PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PetscLimiterSetUp - Construct data structures for the `PetscLimiter`
206: Collective
208: Input Parameter:
209: . lim - the `PetscLimiter` object to setup
211: Level: intermediate
213: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscLimiterDestroy()`
214: @*/
215: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
216: {
217: PetscFunctionBegin;
219: PetscTryTypeMethod(lim, setup);
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: PetscLimiterDestroy - Destroys a `PetscLimiter` object
226: Collective
228: Input Parameter:
229: . lim - the `PetscLimiter` object to destroy
231: Level: beginner
233: .seealso: `PetscLimiter`, `PetscLimiterView()`
234: @*/
235: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
236: {
237: PetscFunctionBegin;
238: if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);
241: if (--((PetscObject)*lim)->refct > 0) {
242: *lim = NULL;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
245: ((PetscObject)*lim)->refct = 0;
247: PetscTryTypeMethod(*lim, destroy);
248: PetscCall(PetscHeaderDestroy(lim));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@
253: PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.
255: Collective
257: Input Parameter:
258: . comm - The communicator for the `PetscLimiter` object
260: Output Parameter:
261: . lim - The `PetscLimiter` object
263: Level: beginner
265: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
266: @*/
267: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
268: {
269: PetscLimiter l;
271: PetscFunctionBegin;
272: PetscAssertPointer(lim, 2);
273: PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
274: PetscCall(PetscFVInitializePackage());
276: PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));
278: *lim = l;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: PetscLimiterLimit - Limit the flux
285: Input Parameters:
286: + lim - The `PetscLimiter`
287: - flim - The input field
289: Output Parameter:
290: . phi - The limited field
292: Level: beginner
294: Note:
295: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
296: .vb
297: The classical flux-limited formulation is psi(r) where
299: r = (u[0] - u[-1]) / (u[1] - u[0])
301: The second order TVD region is bounded by
303: psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
305: where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
306: phi(r)(r+1)/2 in which the reconstructed interface values are
308: u(v) = u[0] + phi(r) (grad u)[0] v
310: where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
312: 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))
314: For a nicer symmetric formulation, rewrite in terms of
316: f = (u[0] - u[-1]) / (u[1] - u[-1])
318: where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
320: phi(r) = phi(1/r)
322: becomes
324: w(f) = w(1-f).
326: The limiters below implement this final form w(f). The reference methods are
328: w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
329: .ve
331: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
332: @*/
333: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
334: {
335: PetscFunctionBegin;
337: PetscAssertPointer(phi, 3);
338: PetscUseTypeMethod(lim, limit, flim, phi);
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
343: {
344: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
346: PetscFunctionBegin;
347: PetscCall(PetscFree(l));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
352: {
353: PetscViewerFormat format;
355: PetscFunctionBegin;
356: PetscCall(PetscViewerGetFormat(viewer, &format));
357: PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
362: {
363: PetscBool isascii;
365: PetscFunctionBegin;
368: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
369: if (isascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
374: {
375: PetscFunctionBegin;
376: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
381: {
382: PetscFunctionBegin;
383: lim->ops->view = PetscLimiterView_Sin;
384: lim->ops->destroy = PetscLimiterDestroy_Sin;
385: lim->ops->limit = PetscLimiterLimit_Sin;
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*MC
390: PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation
392: Level: intermediate
394: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
395: M*/
397: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
398: {
399: PetscLimiter_Sin *l;
401: PetscFunctionBegin;
403: PetscCall(PetscNew(&l));
404: lim->data = l;
406: PetscCall(PetscLimiterInitialize_Sin(lim));
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
411: {
412: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
414: PetscFunctionBegin;
415: PetscCall(PetscFree(l));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
420: {
421: PetscViewerFormat format;
423: PetscFunctionBegin;
424: PetscCall(PetscViewerGetFormat(viewer, &format));
425: PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
430: {
431: PetscBool isascii;
433: PetscFunctionBegin;
436: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
437: if (isascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
442: {
443: PetscFunctionBegin;
444: *phi = 0.0;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
449: {
450: PetscFunctionBegin;
451: lim->ops->view = PetscLimiterView_Zero;
452: lim->ops->destroy = PetscLimiterDestroy_Zero;
453: lim->ops->limit = PetscLimiterLimit_Zero;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*MC
458: PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation
460: Level: intermediate
462: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
463: M*/
465: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
466: {
467: PetscLimiter_Zero *l;
469: PetscFunctionBegin;
471: PetscCall(PetscNew(&l));
472: lim->data = l;
474: PetscCall(PetscLimiterInitialize_Zero(lim));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
479: {
480: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
482: PetscFunctionBegin;
483: PetscCall(PetscFree(l));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
488: {
489: PetscViewerFormat format;
491: PetscFunctionBegin;
492: PetscCall(PetscViewerGetFormat(viewer, &format));
493: PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
498: {
499: PetscBool isascii;
501: PetscFunctionBegin;
504: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
505: if (isascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
510: {
511: PetscFunctionBegin;
512: *phi = 1.0;
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
517: {
518: PetscFunctionBegin;
519: lim->ops->view = PetscLimiterView_None;
520: lim->ops->destroy = PetscLimiterDestroy_None;
521: lim->ops->limit = PetscLimiterLimit_None;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*MC
526: PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation
528: Level: intermediate
530: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
531: M*/
533: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
534: {
535: PetscLimiter_None *l;
537: PetscFunctionBegin;
539: PetscCall(PetscNew(&l));
540: lim->data = l;
542: PetscCall(PetscLimiterInitialize_None(lim));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
547: {
548: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
550: PetscFunctionBegin;
551: PetscCall(PetscFree(l));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
556: {
557: PetscViewerFormat format;
559: PetscFunctionBegin;
560: PetscCall(PetscViewerGetFormat(viewer, &format));
561: PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
566: {
567: PetscBool isascii;
569: PetscFunctionBegin;
572: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
573: if (isascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
578: {
579: PetscFunctionBegin;
580: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
585: {
586: PetscFunctionBegin;
587: lim->ops->view = PetscLimiterView_Minmod;
588: lim->ops->destroy = PetscLimiterDestroy_Minmod;
589: lim->ops->limit = PetscLimiterLimit_Minmod;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*MC
594: PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation
596: Level: intermediate
598: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
599: M*/
601: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
602: {
603: PetscLimiter_Minmod *l;
605: PetscFunctionBegin;
607: PetscCall(PetscNew(&l));
608: lim->data = l;
610: PetscCall(PetscLimiterInitialize_Minmod(lim));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
615: {
616: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
618: PetscFunctionBegin;
619: PetscCall(PetscFree(l));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
624: {
625: PetscViewerFormat format;
627: PetscFunctionBegin;
628: PetscCall(PetscViewerGetFormat(viewer, &format));
629: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
634: {
635: PetscBool isascii;
637: PetscFunctionBegin;
640: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
641: if (isascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
646: {
647: PetscFunctionBegin;
648: *phi = PetscMax(0, 4 * f * (1 - f));
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
653: {
654: PetscFunctionBegin;
655: lim->ops->view = PetscLimiterView_VanLeer;
656: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
657: lim->ops->limit = PetscLimiterLimit_VanLeer;
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*MC
662: PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation
664: Level: intermediate
666: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
667: M*/
669: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
670: {
671: PetscLimiter_VanLeer *l;
673: PetscFunctionBegin;
675: PetscCall(PetscNew(&l));
676: lim->data = l;
678: PetscCall(PetscLimiterInitialize_VanLeer(lim));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
683: {
684: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
686: PetscFunctionBegin;
687: PetscCall(PetscFree(l));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
692: {
693: PetscViewerFormat format;
695: PetscFunctionBegin;
696: PetscCall(PetscViewerGetFormat(viewer, &format));
697: PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
702: {
703: PetscBool isascii;
705: PetscFunctionBegin;
708: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
709: if (isascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
714: {
715: PetscFunctionBegin;
716: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
721: {
722: PetscFunctionBegin;
723: lim->ops->view = PetscLimiterView_VanAlbada;
724: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
725: lim->ops->limit = PetscLimiterLimit_VanAlbada;
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*MC
730: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation
732: Level: intermediate
734: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
735: M*/
737: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
738: {
739: PetscLimiter_VanAlbada *l;
741: PetscFunctionBegin;
743: PetscCall(PetscNew(&l));
744: lim->data = l;
746: PetscCall(PetscLimiterInitialize_VanAlbada(lim));
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
751: {
752: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
754: PetscFunctionBegin;
755: PetscCall(PetscFree(l));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
760: {
761: PetscViewerFormat format;
763: PetscFunctionBegin;
764: PetscCall(PetscViewerGetFormat(viewer, &format));
765: PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
770: {
771: PetscBool isascii;
773: PetscFunctionBegin;
776: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
777: if (isascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
782: {
783: PetscFunctionBegin;
784: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
789: {
790: PetscFunctionBegin;
791: lim->ops->view = PetscLimiterView_Superbee;
792: lim->ops->destroy = PetscLimiterDestroy_Superbee;
793: lim->ops->limit = PetscLimiterLimit_Superbee;
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*MC
798: PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation
800: Level: intermediate
802: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
803: M*/
805: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
806: {
807: PetscLimiter_Superbee *l;
809: PetscFunctionBegin;
811: PetscCall(PetscNew(&l));
812: lim->data = l;
814: PetscCall(PetscLimiterInitialize_Superbee(lim));
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
819: {
820: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
822: PetscFunctionBegin;
823: PetscCall(PetscFree(l));
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
828: {
829: PetscViewerFormat format;
831: PetscFunctionBegin;
832: PetscCall(PetscViewerGetFormat(viewer, &format));
833: PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
838: {
839: PetscBool isascii;
841: PetscFunctionBegin;
844: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
845: if (isascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /* aka Barth-Jespersen */
850: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
851: {
852: PetscFunctionBegin;
853: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
858: {
859: PetscFunctionBegin;
860: lim->ops->view = PetscLimiterView_MC;
861: lim->ops->destroy = PetscLimiterDestroy_MC;
862: lim->ops->limit = PetscLimiterLimit_MC;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*MC
867: PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation
869: Level: intermediate
871: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
872: M*/
874: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
875: {
876: PetscLimiter_MC *l;
878: PetscFunctionBegin;
880: PetscCall(PetscNew(&l));
881: lim->data = l;
883: PetscCall(PetscLimiterInitialize_MC(lim));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: PetscClassId PETSCFV_CLASSID = 0;
889: PetscFunctionList PetscFVList = NULL;
890: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
892: /*@C
893: PetscFVRegister - Adds a new `PetscFV` implementation
895: Not Collective, No Fortran Support
897: Input Parameters:
898: + sname - The name of a new user-defined creation routine
899: - function - The creation routine itself
901: Example Usage:
902: .vb
903: PetscFVRegister("my_fv", MyPetscFVCreate);
904: .ve
906: Then, your PetscFV type can be chosen with the procedural interface via
907: .vb
908: PetscFVCreate(MPI_Comm, PetscFV *);
909: PetscFVSetType(PetscFV, "my_fv");
910: .ve
911: or at runtime via the option
912: .vb
913: -petscfv_type my_fv
914: .ve
916: Level: advanced
918: Note:
919: `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs
921: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
922: @*/
923: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
924: {
925: PetscFunctionBegin;
926: PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: /*@
931: PetscFVSetType - Builds a particular `PetscFV`
933: Collective
935: Input Parameters:
936: + fvm - The `PetscFV` object
937: - name - The type of FVM space
939: Options Database Key:
940: . -petscfv_type type - Sets the `PetscFVType`; use -help for a list of available types
942: Level: intermediate
944: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
945: @*/
946: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
947: {
948: PetscErrorCode (*r)(PetscFV);
949: PetscBool match;
951: PetscFunctionBegin;
953: PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
954: if (match) PetscFunctionReturn(PETSC_SUCCESS);
956: PetscCall(PetscFVRegisterAll());
957: PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
958: PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
960: PetscTryTypeMethod(fvm, destroy);
961: fvm->ops->destroy = NULL;
963: PetscCall((*r)(fvm));
964: PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@
969: PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.
971: Not Collective
973: Input Parameter:
974: . fvm - The `PetscFV`
976: Output Parameter:
977: . name - The `PetscFVType` name
979: Level: intermediate
981: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
982: @*/
983: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
984: {
985: PetscFunctionBegin;
987: PetscAssertPointer(name, 2);
988: PetscCall(PetscFVRegisterAll());
989: *name = ((PetscObject)fvm)->type_name;
990: PetscFunctionReturn(PETSC_SUCCESS);
991: }
993: /*@
994: PetscFVViewFromOptions - View a `PetscFV` based on values in the options database
996: Collective
998: Input Parameters:
999: + A - the `PetscFV` object
1000: . obj - Optional object that provides the options prefix
1001: - name - command line option name
1003: Level: intermediate
1005: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
1006: @*/
1007: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1008: {
1009: PetscFunctionBegin;
1011: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: /*@
1016: PetscFVView - Views a `PetscFV`
1018: Collective
1020: Input Parameters:
1021: + fvm - the `PetscFV` object to view
1022: - v - the viewer
1024: Level: beginner
1026: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1027: @*/
1028: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1029: {
1030: PetscFunctionBegin;
1032: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1033: PetscTryTypeMethod(fvm, view, v);
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: /*@
1038: PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database
1040: Collective
1042: Input Parameter:
1043: . fvm - the `PetscFV` object to set options for
1045: Options Database Key:
1046: . -petscfv_compute_gradients (true|false) - Determines whether cell gradients are calculated
1048: Level: intermediate
1050: .seealso: `PetscFV`, `PetscFVView()`
1051: @*/
1052: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1053: {
1054: const char *defaultType;
1055: char name[256];
1056: PetscBool flg;
1058: PetscFunctionBegin;
1060: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1061: else defaultType = ((PetscObject)fvm)->type_name;
1062: PetscCall(PetscFVRegisterAll());
1064: PetscObjectOptionsBegin((PetscObject)fvm);
1065: PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg));
1066: if (flg) PetscCall(PetscFVSetType(fvm, name));
1067: else if (!((PetscObject)fvm)->type_name) PetscCall(PetscFVSetType(fvm, defaultType));
1068: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1069: PetscTryTypeMethod(fvm, setfromoptions);
1070: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1071: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1072: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1073: PetscOptionsEnd();
1074: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@
1079: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1081: Collective
1083: Input Parameter:
1084: . fvm - the `PetscFV` object to setup
1086: Level: intermediate
1088: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1089: @*/
1090: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1091: {
1092: PetscFunctionBegin;
1094: PetscCall(PetscLimiterSetUp(fvm->limiter));
1095: PetscTryTypeMethod(fvm, setup);
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: /*@
1100: PetscFVDestroy - Destroys a `PetscFV` object
1102: Collective
1104: Input Parameter:
1105: . fvm - the `PetscFV` object to destroy
1107: Level: beginner
1109: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1110: @*/
1111: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1112: {
1113: PetscInt i;
1115: PetscFunctionBegin;
1116: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1119: if (--((PetscObject)*fvm)->refct > 0) {
1120: *fvm = NULL;
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1123: ((PetscObject)*fvm)->refct = 0;
1125: for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1126: PetscCall(PetscFree((*fvm)->componentNames));
1127: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1128: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1129: PetscCall(PetscFree((*fvm)->fluxWork));
1130: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1131: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1133: PetscTryTypeMethod(*fvm, destroy);
1134: PetscCall(PetscHeaderDestroy(fvm));
1135: PetscFunctionReturn(PETSC_SUCCESS);
1136: }
1138: /*@
1139: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1141: Collective
1143: Input Parameter:
1144: . comm - The communicator for the `PetscFV` object
1146: Output Parameter:
1147: . fvm - The `PetscFV` object
1149: Level: beginner
1151: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1152: @*/
1153: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1154: {
1155: PetscFV f;
1157: PetscFunctionBegin;
1158: PetscAssertPointer(fvm, 2);
1159: PetscCall(PetscFVInitializePackage());
1161: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1162: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1163: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1164: f->numComponents = 1;
1165: f->dim = 0;
1166: f->computeGradients = PETSC_FALSE;
1167: f->fluxWork = NULL;
1168: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1170: *fvm = f;
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /*@
1175: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1177: Logically Collective
1179: Input Parameters:
1180: + fvm - the `PetscFV` object
1181: - lim - The `PetscLimiter`
1183: Level: intermediate
1185: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1186: @*/
1187: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1188: {
1189: PetscFunctionBegin;
1192: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1193: PetscCall(PetscObjectReference((PetscObject)lim));
1194: fvm->limiter = lim;
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*@
1199: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1201: Not Collective
1203: Input Parameter:
1204: . fvm - the `PetscFV` object
1206: Output Parameter:
1207: . lim - The `PetscLimiter`
1209: Level: intermediate
1211: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1212: @*/
1213: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1214: {
1215: PetscFunctionBegin;
1217: PetscAssertPointer(lim, 2);
1218: *lim = fvm->limiter;
1219: PetscFunctionReturn(PETSC_SUCCESS);
1220: }
1222: /*@
1223: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1225: Logically Collective
1227: Input Parameters:
1228: + fvm - the `PetscFV` object
1229: - comp - The number of components
1231: Level: intermediate
1233: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1234: @*/
1235: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1236: {
1237: PetscFunctionBegin;
1239: if (fvm->numComponents != comp) {
1240: PetscInt i;
1242: for (i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1243: PetscCall(PetscFree(fvm->componentNames));
1244: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1245: }
1246: fvm->numComponents = comp;
1247: PetscCall(PetscFree(fvm->fluxWork));
1248: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }
1252: /*@
1253: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1255: Not Collective
1257: Input Parameter:
1258: . fvm - the `PetscFV` object
1260: Output Parameter:
1261: . comp - The number of components
1263: Level: intermediate
1265: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1266: @*/
1267: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1268: {
1269: PetscFunctionBegin;
1271: PetscAssertPointer(comp, 2);
1272: *comp = fvm->numComponents;
1273: PetscFunctionReturn(PETSC_SUCCESS);
1274: }
1276: /*@
1277: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1279: Logically Collective
1281: Input Parameters:
1282: + fvm - the `PetscFV` object
1283: . comp - the component number
1284: - name - the component name
1286: Level: intermediate
1288: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1289: @*/
1290: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1291: {
1292: PetscFunctionBegin;
1293: PetscCall(PetscFree(fvm->componentNames[comp]));
1294: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1295: PetscFunctionReturn(PETSC_SUCCESS);
1296: }
1298: /*@
1299: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1301: Logically Collective
1303: Input Parameters:
1304: + fvm - the `PetscFV` object
1305: - comp - the component number
1307: Output Parameter:
1308: . name - the component name
1310: Level: intermediate
1312: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1313: @*/
1314: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char *name[])
1315: {
1316: PetscFunctionBegin;
1317: *name = fvm->componentNames[comp];
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: /*@
1322: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1324: Logically Collective
1326: Input Parameters:
1327: + fvm - the `PetscFV` object
1328: - dim - The spatial dimension
1330: Level: intermediate
1332: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
1333: @*/
1334: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1335: {
1336: PetscFunctionBegin;
1338: fvm->dim = dim;
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: /*@
1343: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1345: Not Collective
1347: Input Parameter:
1348: . fvm - the `PetscFV` object
1350: Output Parameter:
1351: . dim - The spatial dimension
1353: Level: intermediate
1355: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1356: @*/
1357: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1358: {
1359: PetscFunctionBegin;
1361: PetscAssertPointer(dim, 2);
1362: *dim = fvm->dim;
1363: PetscFunctionReturn(PETSC_SUCCESS);
1364: }
1366: /*@
1367: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1369: Logically Collective
1371: Input Parameters:
1372: + fvm - the `PetscFV` object
1373: - computeGradients - Flag to compute cell gradients
1375: Level: intermediate
1377: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1378: @*/
1379: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1380: {
1381: PetscFunctionBegin;
1383: fvm->computeGradients = computeGradients;
1384: PetscFunctionReturn(PETSC_SUCCESS);
1385: }
1387: /*@
1388: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1390: Not Collective
1392: Input Parameter:
1393: . fvm - the `PetscFV` object
1395: Output Parameter:
1396: . computeGradients - Flag to compute cell gradients
1398: Level: intermediate
1400: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1401: @*/
1402: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1403: {
1404: PetscFunctionBegin;
1406: PetscAssertPointer(computeGradients, 2);
1407: *computeGradients = fvm->computeGradients;
1408: PetscFunctionReturn(PETSC_SUCCESS);
1409: }
1411: /*@
1412: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1414: Logically Collective
1416: Input Parameters:
1417: + fvm - the `PetscFV` object
1418: - q - The `PetscQuadrature`
1420: Level: intermediate
1422: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1423: @*/
1424: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1425: {
1426: PetscFunctionBegin;
1428: PetscCall(PetscObjectReference((PetscObject)q));
1429: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1430: fvm->quadrature = q;
1431: PetscFunctionReturn(PETSC_SUCCESS);
1432: }
1434: /*@
1435: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1437: Not Collective
1439: Input Parameter:
1440: . fvm - the `PetscFV` object
1442: Output Parameter:
1443: . q - The `PetscQuadrature`
1445: Level: intermediate
1447: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1448: @*/
1449: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1450: {
1451: PetscFunctionBegin;
1453: PetscAssertPointer(q, 2);
1454: if (!fvm->quadrature) {
1455: /* Create default 1-point quadrature */
1456: PetscReal *points, *weights;
1458: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1459: PetscCall(PetscCalloc1(fvm->dim, &points));
1460: PetscCall(PetscMalloc1(1, &weights));
1461: weights[0] = 1.0;
1462: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1463: }
1464: *q = fvm->quadrature;
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: /*@
1469: PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`
1471: Not Collective
1473: Input Parameters:
1474: + fvm - The `PetscFV` object
1475: - ct - The `DMPolytopeType` for the cell
1477: Level: intermediate
1479: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1480: @*/
1481: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1482: {
1483: DM K;
1484: PetscInt dim, Nc;
1486: PetscFunctionBegin;
1487: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1488: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1489: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1490: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1491: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1492: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1493: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1494: PetscCall(DMDestroy(&K));
1495: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1496: // Should we be using PetscFVGetQuadrature() here?
1497: for (PetscInt c = 0; c < Nc; ++c) {
1498: PetscQuadrature qc;
1499: PetscReal *points, *weights;
1501: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1502: PetscCall(PetscCalloc1(dim, &points));
1503: PetscCall(PetscCalloc1(Nc, &weights));
1504: weights[c] = 1.0;
1505: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1506: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1507: PetscCall(PetscQuadratureDestroy(&qc));
1508: }
1509: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: /*@
1514: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1516: Not Collective
1518: Input Parameter:
1519: . fvm - The `PetscFV` object
1521: Output Parameter:
1522: . sp - The `PetscDualSpace` object
1524: Level: intermediate
1526: Developer Notes:
1527: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1529: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1530: @*/
1531: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1532: {
1533: PetscFunctionBegin;
1535: PetscAssertPointer(sp, 2);
1536: if (!fvm->dualSpace) {
1537: PetscInt dim;
1539: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1540: PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1541: }
1542: *sp = fvm->dualSpace;
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: /*@
1547: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1549: Not Collective
1551: Input Parameters:
1552: + fvm - The `PetscFV` object
1553: - sp - The `PetscDualSpace` object
1555: Level: intermediate
1557: Note:
1558: A simple dual space is provided automatically, and the user typically will not need to override it.
1560: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1561: @*/
1562: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1563: {
1564: PetscFunctionBegin;
1567: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1568: fvm->dualSpace = sp;
1569: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1570: PetscFunctionReturn(PETSC_SUCCESS);
1571: }
1573: /*@C
1574: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1576: Not Collective
1578: Input Parameter:
1579: . fvm - The `PetscFV` object
1581: Output Parameter:
1582: . T - The basis function values and derivatives at quadrature points
1584: Level: intermediate
1586: Note:
1587: .vb
1588: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1589: 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
1590: 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
1591: .ve
1593: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1594: @*/
1595: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1596: {
1597: PetscInt npoints;
1598: const PetscReal *points;
1600: PetscFunctionBegin;
1602: PetscAssertPointer(T, 2);
1603: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1604: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1605: *T = fvm->T;
1606: PetscFunctionReturn(PETSC_SUCCESS);
1607: }
1609: /*@C
1610: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1612: Not Collective
1614: Input Parameters:
1615: + fvm - The `PetscFV` object
1616: . nrepl - The number of replicas
1617: . npoints - The number of tabulation points in a replica
1618: . points - The tabulation point coordinates
1619: - K - The order of derivative to tabulate
1621: Output Parameter:
1622: . T - The basis function values and derivative at tabulation points
1624: Level: intermediate
1626: Note:
1627: .vb
1628: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1629: 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
1630: 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
1631: .ve
1633: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1634: @*/
1635: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1636: {
1637: PetscInt pdim; // Dimension of approximation space P
1638: PetscInt cdim; // Spatial dimension
1639: PetscInt Nc; // Field components
1640: PetscInt k, p, d, c, e;
1642: PetscFunctionBegin;
1643: if (!npoints || K < 0) {
1644: *T = NULL;
1645: PetscFunctionReturn(PETSC_SUCCESS);
1646: }
1648: PetscAssertPointer(points, 4);
1649: PetscAssertPointer(T, 6);
1650: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1651: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1652: pdim = Nc;
1653: PetscCall(PetscMalloc1(1, T));
1654: (*T)->K = !cdim ? 0 : K;
1655: (*T)->Nr = nrepl;
1656: (*T)->Np = npoints;
1657: (*T)->Nb = pdim;
1658: (*T)->Nc = Nc;
1659: (*T)->cdim = cdim;
1660: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1661: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1662: if (K >= 0) {
1663: for (p = 0; p < nrepl * npoints; ++p)
1664: for (d = 0; d < pdim; ++d)
1665: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1666: }
1667: if (K >= 1) {
1668: for (p = 0; p < nrepl * npoints; ++p)
1669: for (d = 0; d < pdim; ++d)
1670: for (c = 0; c < Nc; ++c)
1671: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1672: }
1673: if (K >= 2) {
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 * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1678: }
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: /*@
1683: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1685: Input Parameters:
1686: + fvm - The `PetscFV` object
1687: . numFaces - The number of cell faces which are not constrained
1688: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1690: Output Parameter:
1691: . grad - the gradient
1693: Level: advanced
1695: .seealso: `PetscFV`, `PetscFVCreate()`
1696: @*/
1697: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1698: {
1699: PetscFunctionBegin;
1701: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: /*@C
1706: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1708: Not Collective
1710: Input Parameters:
1711: + fvm - The `PetscFV` object for the field being integrated
1712: . prob - The `PetscDS` specifying the discretizations and continuum functions
1713: . field - The field being integrated
1714: . Nf - The number of faces in the chunk
1715: . fgeom - The face geometry for each face in the chunk
1716: . neighborVol - The volume for each pair of cells in the chunk
1717: . uL - The state from the cell on the left
1718: - uR - The state from the cell on the right
1720: Output Parameters:
1721: + fluxL - the left fluxes for each face
1722: - fluxR - the right fluxes for each face
1724: Level: developer
1726: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1727: @*/
1728: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1729: {
1730: PetscFunctionBegin;
1732: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1733: PetscFunctionReturn(PETSC_SUCCESS);
1734: }
1736: /*@
1737: PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.
1739: Input Parameter:
1740: . fv - The initial `PetscFV`
1742: Output Parameter:
1743: . fvNew - A clone of the `PetscFV`
1745: Level: advanced
1747: Notes:
1748: This is typically used to change the number of components.
1750: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1751: @*/
1752: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1753: {
1754: PetscDualSpace Q;
1755: DM K;
1756: PetscQuadrature q;
1757: PetscInt Nc, cdim;
1759: PetscFunctionBegin;
1760: PetscCall(PetscFVGetDualSpace(fv, &Q));
1761: PetscCall(PetscFVGetQuadrature(fv, &q));
1762: PetscCall(PetscDualSpaceGetDM(Q, &K));
1764: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1765: PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1766: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1767: PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1768: PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1769: PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1770: PetscCall(PetscFVSetQuadrature(*fvNew, q));
1771: PetscFunctionReturn(PETSC_SUCCESS);
1772: }
1774: /*@
1775: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1776: smaller copies.
1778: Input Parameter:
1779: . fv - The initial `PetscFV`
1781: Output Parameter:
1782: . fvRef - The refined `PetscFV`
1784: Level: advanced
1786: Notes:
1787: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1788: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1789: interpolation between regularly refined meshes.
1791: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1792: @*/
1793: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1794: {
1795: PetscDualSpace Q, Qref;
1796: DM K, Kref;
1797: PetscQuadrature q, qref;
1798: DMPolytopeType ct;
1799: DMPlexTransform tr;
1800: PetscReal *v0;
1801: PetscReal *jac, *invjac;
1802: PetscInt numComp, numSubelements, s;
1804: PetscFunctionBegin;
1805: PetscCall(PetscFVGetDualSpace(fv, &Q));
1806: PetscCall(PetscFVGetQuadrature(fv, &q));
1807: PetscCall(PetscDualSpaceGetDM(Q, &K));
1808: /* Create dual space */
1809: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1810: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1811: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1812: PetscCall(DMDestroy(&Kref));
1813: PetscCall(PetscDualSpaceSetUp(Qref));
1814: /* Create volume */
1815: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1816: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1817: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1818: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1819: PetscCall(PetscFVSetUp(*fvRef));
1820: /* Create quadrature */
1821: PetscCall(DMPlexGetCellType(K, 0, &ct));
1822: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1823: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1824: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1825: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1826: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1827: for (s = 0; s < numSubelements; ++s) {
1828: PetscQuadrature qs;
1829: const PetscReal *points, *weights;
1830: PetscReal *p, *w;
1831: PetscInt dim, Nc, npoints, np;
1833: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1834: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1835: np = npoints / numSubelements;
1836: PetscCall(PetscMalloc1(np * dim, &p));
1837: PetscCall(PetscMalloc1(np * Nc, &w));
1838: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1839: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1840: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1841: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1842: PetscCall(PetscQuadratureDestroy(&qs));
1843: }
1844: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1845: PetscCall(DMPlexTransformDestroy(&tr));
1846: PetscCall(PetscQuadratureDestroy(&qref));
1847: PetscCall(PetscDualSpaceDestroy(&Qref));
1848: PetscFunctionReturn(PETSC_SUCCESS);
1849: }
1851: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1852: {
1853: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1855: PetscFunctionBegin;
1856: PetscCall(PetscFree(b));
1857: PetscFunctionReturn(PETSC_SUCCESS);
1858: }
1860: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1861: {
1862: PetscInt Nc, c;
1863: PetscViewerFormat format;
1865: PetscFunctionBegin;
1866: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1867: PetscCall(PetscViewerGetFormat(viewer, &format));
1868: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1869: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1870: for (c = 0; c < Nc; c++) {
1871: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1872: }
1873: PetscFunctionReturn(PETSC_SUCCESS);
1874: }
1876: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1877: {
1878: PetscBool isascii;
1880: PetscFunctionBegin;
1883: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1884: if (isascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1885: PetscFunctionReturn(PETSC_SUCCESS);
1886: }
1888: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1889: {
1890: PetscInt dim;
1892: PetscFunctionBegin;
1893: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1894: for (PetscInt f = 0; f < numFaces; ++f) {
1895: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1896: }
1897: PetscFunctionReturn(PETSC_SUCCESS);
1898: }
1900: /*
1901: neighborVol[f*2+0] contains the left geom
1902: neighborVol[f*2+1] contains the right geom
1903: */
1904: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1905: {
1906: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1907: void *rctx;
1908: PetscScalar *flux = fvm->fluxWork;
1909: const PetscScalar *constants;
1910: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1912: PetscFunctionBegin;
1913: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1914: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1915: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1916: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1917: PetscCall(PetscDSGetContext(prob, field, &rctx));
1918: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1919: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1920: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1921: for (f = 0; f < Nf; ++f) {
1922: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1923: for (d = 0; d < pdim; ++d) {
1924: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1925: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1926: }
1927: }
1928: PetscFunctionReturn(PETSC_SUCCESS);
1929: }
1931: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1932: {
1933: PetscFunctionBegin;
1934: fvm->ops->setfromoptions = NULL;
1935: fvm->ops->view = PetscFVView_Upwind;
1936: fvm->ops->destroy = PetscFVDestroy_Upwind;
1937: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1938: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1939: PetscFunctionReturn(PETSC_SUCCESS);
1940: }
1942: /*MC
1943: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1945: Level: intermediate
1947: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1948: M*/
1950: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1951: {
1952: PetscFV_Upwind *b;
1954: PetscFunctionBegin;
1956: PetscCall(PetscNew(&b));
1957: fvm->data = b;
1959: PetscCall(PetscFVInitialize_Upwind(fvm));
1960: PetscFunctionReturn(PETSC_SUCCESS);
1961: }
1963: #include <petscblaslapack.h>
1965: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1966: {
1967: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1969: PetscFunctionBegin;
1970: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1971: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1972: PetscCall(PetscFree(ls));
1973: PetscFunctionReturn(PETSC_SUCCESS);
1974: }
1976: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1977: {
1978: PetscInt Nc, c;
1979: PetscViewerFormat format;
1981: PetscFunctionBegin;
1982: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1983: PetscCall(PetscViewerGetFormat(viewer, &format));
1984: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1985: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1986: for (c = 0; c < Nc; c++) {
1987: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1988: }
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1993: {
1994: PetscBool isascii;
1996: PetscFunctionBegin;
1999: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2000: if (isascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2001: PetscFunctionReturn(PETSC_SUCCESS);
2002: }
2004: /* Overwrites A. Can only handle full-rank problems with m>=n */
2005: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2006: {
2007: PetscBool debug = PETSC_FALSE;
2008: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2009: PetscScalar *R, *Q, *Aback, Alpha;
2011: PetscFunctionBegin;
2012: if (debug) {
2013: PetscCall(PetscMalloc1(m * n, &Aback));
2014: PetscCall(PetscArraycpy(Aback, A, m * n));
2015: }
2017: PetscCall(PetscBLASIntCast(m, &M));
2018: PetscCall(PetscBLASIntCast(n, &N));
2019: PetscCall(PetscBLASIntCast(mstride, &lda));
2020: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2021: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2022: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2023: PetscCall(PetscFPTrapPop());
2024: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2025: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2027: /* Extract an explicit representation of Q */
2028: Q = Ainv;
2029: PetscCall(PetscArraycpy(Q, A, mstride * n));
2030: K = N; /* full rank */
2031: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2032: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2034: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2035: Alpha = 1.0;
2036: ldb = lda;
2037: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2038: /* Ainv is Q, overwritten with inverse */
2040: if (debug) { /* Check that pseudo-inverse worked */
2041: PetscScalar Beta = 0.0;
2042: PetscBLASInt ldc;
2043: K = N;
2044: ldc = N;
2045: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2046: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2047: PetscCall(PetscFree(Aback));
2048: }
2049: PetscFunctionReturn(PETSC_SUCCESS);
2050: }
2052: /* Overwrites A. Can handle degenerate problems and m<n. */
2053: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2054: {
2055: PetscScalar *Brhs;
2056: PetscScalar *tmpwork;
2057: PetscReal rcond;
2058: #if defined(PETSC_USE_COMPLEX)
2059: PetscInt rworkSize;
2060: PetscReal *rwork, *rtau;
2061: #endif
2062: PetscInt i, j, maxmn;
2063: PetscBLASInt M, N, lda, ldb, ldwork;
2064: PetscBLASInt nrhs, irank, info;
2066: PetscFunctionBegin;
2067: /* initialize to identity */
2068: tmpwork = work;
2069: Brhs = Ainv;
2070: maxmn = PetscMax(m, n);
2071: for (j = 0; j < maxmn; j++) {
2072: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2073: }
2075: PetscCall(PetscBLASIntCast(m, &M));
2076: PetscCall(PetscBLASIntCast(n, &N));
2077: PetscCall(PetscBLASIntCast(mstride, &lda));
2078: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2079: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2080: rcond = -1;
2081: nrhs = M;
2082: #if defined(PETSC_USE_COMPLEX)
2083: rworkSize = 5 * PetscMin(M, N);
2084: PetscCall(PetscMalloc1(rworkSize, &rwork));
2085: PetscCall(PetscMalloc1(PetscMin(M, N), &rtau));
2086: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2087: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, rtau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2088: PetscCall(PetscFPTrapPop());
2089: PetscCall(PetscFree(rwork));
2090: for (i = 0; i < PetscMin(M, N); i++) tau[i] = rtau[i];
2091: PetscCall(PetscFree(rtau));
2092: #else
2093: nrhs = M;
2094: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2095: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, tau, &rcond, &irank, tmpwork, &ldwork, &info));
2096: PetscCall(PetscFPTrapPop());
2097: #endif
2098: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2099: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2100: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two collinear points");
2101: PetscFunctionReturn(PETSC_SUCCESS);
2102: }
2104: #if 0
2105: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2106: {
2107: PetscReal grad[2] = {0, 0};
2108: const PetscInt *faces;
2109: PetscInt numFaces, f;
2111: PetscFunctionBegin;
2112: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2113: PetscCall(DMPlexGetCone(dm, cell, &faces));
2114: for (f = 0; f < numFaces; ++f) {
2115: const PetscInt *fcells;
2116: const CellGeom *cg1;
2117: const FaceGeom *fg;
2119: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2120: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2121: for (i = 0; i < 2; ++i) {
2122: PetscScalar du;
2124: if (fcells[i] == c) continue;
2125: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2126: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2127: grad[0] += fg->grad[!i][0] * du;
2128: grad[1] += fg->grad[!i][1] * du;
2129: }
2130: }
2131: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2132: PetscFunctionReturn(PETSC_SUCCESS);
2133: }
2134: #endif
2136: /*
2137: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2139: Input Parameters:
2140: + fvm - The `PetscFV` object
2141: . numFaces - The number of cell faces which are not constrained
2142: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2144: Level: developer
2146: .seealso: `PetscFV`, `PetscFVCreate()`
2147: */
2148: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2149: {
2150: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2151: const PetscBool useSVD = PETSC_TRUE;
2152: const PetscInt maxFaces = ls->maxFaces;
2153: PetscInt dim, f, d;
2155: PetscFunctionBegin;
2156: if (numFaces > maxFaces) {
2157: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2158: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2159: }
2160: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2161: for (f = 0; f < numFaces; ++f) {
2162: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2163: }
2164: /* Overwrites B with garbage, returns Binv in row-major format */
2165: if (useSVD) {
2166: PetscInt maxmn = PetscMax(numFaces, dim);
2167: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2168: /* Binv shaped in column-major, coldim=maxmn.*/
2169: for (f = 0; f < numFaces; ++f) {
2170: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2171: }
2172: } else {
2173: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2174: /* Binv shaped in row-major, rowdim=maxFaces.*/
2175: for (f = 0; f < numFaces; ++f) {
2176: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2177: }
2178: }
2179: PetscFunctionReturn(PETSC_SUCCESS);
2180: }
2182: /*
2183: neighborVol[f*2+0] contains the left geom
2184: neighborVol[f*2+1] contains the right geom
2185: */
2186: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2187: {
2188: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2189: void *rctx;
2190: PetscScalar *flux = fvm->fluxWork;
2191: const PetscScalar *constants;
2192: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2194: PetscFunctionBegin;
2195: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2196: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2197: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2198: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2199: PetscCall(PetscDSGetContext(prob, field, &rctx));
2200: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2201: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2202: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2203: for (f = 0; f < Nf; ++f) {
2204: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2205: for (d = 0; d < pdim; ++d) {
2206: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2207: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2208: }
2209: }
2210: PetscFunctionReturn(PETSC_SUCCESS);
2211: }
2213: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2214: {
2215: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2216: PetscInt dim, m, n, nrhs, minmn, maxmn;
2218: PetscFunctionBegin;
2220: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2221: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2222: ls->maxFaces = maxFaces;
2223: m = ls->maxFaces;
2224: n = dim;
2225: nrhs = ls->maxFaces;
2226: minmn = PetscMin(m, n);
2227: maxmn = PetscMax(m, n);
2228: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2229: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2230: PetscFunctionReturn(PETSC_SUCCESS);
2231: }
2233: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2234: {
2235: PetscFunctionBegin;
2236: fvm->ops->setfromoptions = NULL;
2237: fvm->ops->view = PetscFVView_LeastSquares;
2238: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2239: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2240: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2241: PetscFunctionReturn(PETSC_SUCCESS);
2242: }
2244: /*MC
2245: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2247: Level: intermediate
2249: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2250: M*/
2252: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2253: {
2254: PetscFV_LeastSquares *ls;
2256: PetscFunctionBegin;
2258: PetscCall(PetscNew(&ls));
2259: fvm->data = ls;
2261: ls->maxFaces = -1;
2262: ls->workSize = -1;
2263: ls->B = NULL;
2264: ls->Binv = NULL;
2265: ls->tau = NULL;
2266: ls->work = NULL;
2268: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2269: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2270: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2271: PetscFunctionReturn(PETSC_SUCCESS);
2272: }
2274: /*@
2275: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2277: Not Collective
2279: Input Parameters:
2280: + fvm - The `PetscFV` object
2281: - maxFaces - The maximum number of cell faces
2283: Level: intermediate
2285: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2286: @*/
2287: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2288: {
2289: PetscFunctionBegin;
2291: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2292: PetscFunctionReturn(PETSC_SUCCESS);
2293: }