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 iascii;
365: PetscFunctionBegin;
368: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
369: if (iascii) 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 iascii;
433: PetscFunctionBegin;
436: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
437: if (iascii) 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 iascii;
501: PetscFunctionBegin;
504: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
505: if (iascii) 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 iascii;
569: PetscFunctionBegin;
572: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
573: if (iascii) 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 iascii;
637: PetscFunctionBegin;
640: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
641: if (iascii) 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 iascii;
705: PetscFunctionBegin;
708: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
709: if (iascii) 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 iascii;
773: PetscFunctionBegin;
776: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
777: if (iascii) 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 iascii;
841: PetscFunctionBegin;
844: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
845: if (iascii) 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 <bool> - 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) {
1067: PetscCall(PetscFVSetType(fvm, name));
1068: } else if (!((PetscObject)fvm)->type_name) {
1069: PetscCall(PetscFVSetType(fvm, defaultType));
1070: }
1071: PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1072: PetscTryTypeMethod(fvm, setfromoptions);
1073: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1074: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1075: PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1076: PetscOptionsEnd();
1077: PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: /*@
1082: PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`
1084: Collective
1086: Input Parameter:
1087: . fvm - the `PetscFV` object to setup
1089: Level: intermediate
1091: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1092: @*/
1093: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1094: {
1095: PetscFunctionBegin;
1097: PetscCall(PetscLimiterSetUp(fvm->limiter));
1098: PetscTryTypeMethod(fvm, setup);
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: /*@
1103: PetscFVDestroy - Destroys a `PetscFV` object
1105: Collective
1107: Input Parameter:
1108: . fvm - the `PetscFV` object to destroy
1110: Level: beginner
1112: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1113: @*/
1114: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1115: {
1116: PetscInt i;
1118: PetscFunctionBegin;
1119: if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);
1122: if (--((PetscObject)*fvm)->refct > 0) {
1123: *fvm = NULL;
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1126: ((PetscObject)*fvm)->refct = 0;
1128: for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1129: PetscCall(PetscFree((*fvm)->componentNames));
1130: PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1131: PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1132: PetscCall(PetscFree((*fvm)->fluxWork));
1133: PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1134: PetscCall(PetscTabulationDestroy(&(*fvm)->T));
1136: PetscTryTypeMethod(*fvm, destroy);
1137: PetscCall(PetscHeaderDestroy(fvm));
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: /*@
1142: PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.
1144: Collective
1146: Input Parameter:
1147: . comm - The communicator for the `PetscFV` object
1149: Output Parameter:
1150: . fvm - The `PetscFV` object
1152: Level: beginner
1154: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1155: @*/
1156: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1157: {
1158: PetscFV f;
1160: PetscFunctionBegin;
1161: PetscAssertPointer(fvm, 2);
1162: PetscCall(PetscFVInitializePackage());
1164: PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1165: PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
1166: PetscCall(PetscLimiterCreate(comm, &f->limiter));
1167: f->numComponents = 1;
1168: f->dim = 0;
1169: f->computeGradients = PETSC_FALSE;
1170: f->fluxWork = NULL;
1171: PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));
1173: *fvm = f;
1174: PetscFunctionReturn(PETSC_SUCCESS);
1175: }
1177: /*@
1178: PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`
1180: Logically Collective
1182: Input Parameters:
1183: + fvm - the `PetscFV` object
1184: - lim - The `PetscLimiter`
1186: Level: intermediate
1188: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1189: @*/
1190: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1191: {
1192: PetscFunctionBegin;
1195: PetscCall(PetscLimiterDestroy(&fvm->limiter));
1196: PetscCall(PetscObjectReference((PetscObject)lim));
1197: fvm->limiter = lim;
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: /*@
1202: PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`
1204: Not Collective
1206: Input Parameter:
1207: . fvm - the `PetscFV` object
1209: Output Parameter:
1210: . lim - The `PetscLimiter`
1212: Level: intermediate
1214: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1215: @*/
1216: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1217: {
1218: PetscFunctionBegin;
1220: PetscAssertPointer(lim, 2);
1221: *lim = fvm->limiter;
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: /*@
1226: PetscFVSetNumComponents - Set the number of field components in a `PetscFV`
1228: Logically Collective
1230: Input Parameters:
1231: + fvm - the `PetscFV` object
1232: - comp - The number of components
1234: Level: intermediate
1236: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1237: @*/
1238: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1239: {
1240: PetscFunctionBegin;
1242: if (fvm->numComponents != comp) {
1243: PetscInt i;
1245: for (i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1246: PetscCall(PetscFree(fvm->componentNames));
1247: PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1248: }
1249: fvm->numComponents = comp;
1250: PetscCall(PetscFree(fvm->fluxWork));
1251: PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: PetscFVGetNumComponents - Get the number of field components in a `PetscFV`
1258: Not Collective
1260: Input Parameter:
1261: . fvm - the `PetscFV` object
1263: Output Parameter:
1264: . comp - The number of components
1266: Level: intermediate
1268: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1269: @*/
1270: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1271: {
1272: PetscFunctionBegin;
1274: PetscAssertPointer(comp, 2);
1275: *comp = fvm->numComponents;
1276: PetscFunctionReturn(PETSC_SUCCESS);
1277: }
1279: /*@
1280: PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`
1282: Logically Collective
1284: Input Parameters:
1285: + fvm - the `PetscFV` object
1286: . comp - the component number
1287: - name - the component name
1289: Level: intermediate
1291: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1292: @*/
1293: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1294: {
1295: PetscFunctionBegin;
1296: PetscCall(PetscFree(fvm->componentNames[comp]));
1297: PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: /*@
1302: PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`
1304: Logically Collective
1306: Input Parameters:
1307: + fvm - the `PetscFV` object
1308: - comp - the component number
1310: Output Parameter:
1311: . name - the component name
1313: Level: intermediate
1315: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1316: @*/
1317: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char *name[])
1318: {
1319: PetscFunctionBegin;
1320: *name = fvm->componentNames[comp];
1321: PetscFunctionReturn(PETSC_SUCCESS);
1322: }
1324: /*@
1325: PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`
1327: Logically Collective
1329: Input Parameters:
1330: + fvm - the `PetscFV` object
1331: - dim - The spatial dimension
1333: Level: intermediate
1335: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
1336: @*/
1337: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1338: {
1339: PetscFunctionBegin;
1341: fvm->dim = dim;
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: /*@
1346: PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`
1348: Not Collective
1350: Input Parameter:
1351: . fvm - the `PetscFV` object
1353: Output Parameter:
1354: . dim - The spatial dimension
1356: Level: intermediate
1358: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1359: @*/
1360: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1361: {
1362: PetscFunctionBegin;
1364: PetscAssertPointer(dim, 2);
1365: *dim = fvm->dim;
1366: PetscFunctionReturn(PETSC_SUCCESS);
1367: }
1369: /*@
1370: PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`
1372: Logically Collective
1374: Input Parameters:
1375: + fvm - the `PetscFV` object
1376: - computeGradients - Flag to compute cell gradients
1378: Level: intermediate
1380: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1381: @*/
1382: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1383: {
1384: PetscFunctionBegin;
1386: fvm->computeGradients = computeGradients;
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: /*@
1391: PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`
1393: Not Collective
1395: Input Parameter:
1396: . fvm - the `PetscFV` object
1398: Output Parameter:
1399: . computeGradients - Flag to compute cell gradients
1401: Level: intermediate
1403: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1404: @*/
1405: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1406: {
1407: PetscFunctionBegin;
1409: PetscAssertPointer(computeGradients, 2);
1410: *computeGradients = fvm->computeGradients;
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /*@
1415: PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`
1417: Logically Collective
1419: Input Parameters:
1420: + fvm - the `PetscFV` object
1421: - q - The `PetscQuadrature`
1423: Level: intermediate
1425: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1426: @*/
1427: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1428: {
1429: PetscFunctionBegin;
1431: PetscCall(PetscObjectReference((PetscObject)q));
1432: PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1433: fvm->quadrature = q;
1434: PetscFunctionReturn(PETSC_SUCCESS);
1435: }
1437: /*@
1438: PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`
1440: Not Collective
1442: Input Parameter:
1443: . fvm - the `PetscFV` object
1445: Output Parameter:
1446: . q - The `PetscQuadrature`
1448: Level: intermediate
1450: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1451: @*/
1452: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1453: {
1454: PetscFunctionBegin;
1456: PetscAssertPointer(q, 2);
1457: if (!fvm->quadrature) {
1458: /* Create default 1-point quadrature */
1459: PetscReal *points, *weights;
1461: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1462: PetscCall(PetscCalloc1(fvm->dim, &points));
1463: PetscCall(PetscMalloc1(1, &weights));
1464: weights[0] = 1.0;
1465: PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1466: }
1467: *q = fvm->quadrature;
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: /*@
1472: PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`
1474: Not Collective
1476: Input Parameters:
1477: + fvm - The `PetscFV` object
1478: - ct - The `DMPolytopeType` for the cell
1480: Level: intermediate
1482: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1483: @*/
1484: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1485: {
1486: DM K;
1487: PetscInt dim, Nc;
1489: PetscFunctionBegin;
1490: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1491: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1492: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1493: PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1494: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1495: PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1496: PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1497: PetscCall(DMDestroy(&K));
1498: PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1499: // Should we be using PetscFVGetQuadrature() here?
1500: for (PetscInt c = 0; c < Nc; ++c) {
1501: PetscQuadrature qc;
1502: PetscReal *points, *weights;
1504: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1505: PetscCall(PetscCalloc1(dim, &points));
1506: PetscCall(PetscCalloc1(Nc, &weights));
1507: weights[c] = 1.0;
1508: PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1509: PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1510: PetscCall(PetscQuadratureDestroy(&qc));
1511: }
1512: PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1513: PetscFunctionReturn(PETSC_SUCCESS);
1514: }
1516: /*@
1517: PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`
1519: Not Collective
1521: Input Parameter:
1522: . fvm - The `PetscFV` object
1524: Output Parameter:
1525: . sp - The `PetscDualSpace` object
1527: Level: intermediate
1529: Developer Notes:
1530: There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class
1532: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1533: @*/
1534: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1535: {
1536: PetscFunctionBegin;
1538: PetscAssertPointer(sp, 2);
1539: if (!fvm->dualSpace) {
1540: PetscInt dim;
1542: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1543: PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1544: }
1545: *sp = fvm->dualSpace;
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: /*@
1550: PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product
1552: Not Collective
1554: Input Parameters:
1555: + fvm - The `PetscFV` object
1556: - sp - The `PetscDualSpace` object
1558: Level: intermediate
1560: Note:
1561: A simple dual space is provided automatically, and the user typically will not need to override it.
1563: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1564: @*/
1565: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1566: {
1567: PetscFunctionBegin;
1570: PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1571: fvm->dualSpace = sp;
1572: PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1573: PetscFunctionReturn(PETSC_SUCCESS);
1574: }
1576: /*@C
1577: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1579: Not Collective
1581: Input Parameter:
1582: . fvm - The `PetscFV` object
1584: Output Parameter:
1585: . T - The basis function values and derivatives at quadrature points
1587: Level: intermediate
1589: Note:
1590: .vb
1591: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1592: 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
1593: 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
1594: .ve
1596: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1597: @*/
1598: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1599: {
1600: PetscInt npoints;
1601: const PetscReal *points;
1603: PetscFunctionBegin;
1605: PetscAssertPointer(T, 2);
1606: PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1607: if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1608: *T = fvm->T;
1609: PetscFunctionReturn(PETSC_SUCCESS);
1610: }
1612: /*@C
1613: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1615: Not Collective
1617: Input Parameters:
1618: + fvm - The `PetscFV` object
1619: . nrepl - The number of replicas
1620: . npoints - The number of tabulation points in a replica
1621: . points - The tabulation point coordinates
1622: - K - The order of derivative to tabulate
1624: Output Parameter:
1625: . T - The basis function values and derivative at tabulation points
1627: Level: intermediate
1629: Note:
1630: .vb
1631: T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1632: 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
1633: 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
1634: .ve
1636: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1637: @*/
1638: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1639: {
1640: PetscInt pdim; // Dimension of approximation space P
1641: PetscInt cdim; // Spatial dimension
1642: PetscInt Nc; // Field components
1643: PetscInt k, p, d, c, e;
1645: PetscFunctionBegin;
1646: if (!npoints || K < 0) {
1647: *T = NULL;
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: PetscAssertPointer(points, 4);
1652: PetscAssertPointer(T, 6);
1653: PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1654: PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1655: pdim = Nc;
1656: PetscCall(PetscMalloc1(1, T));
1657: (*T)->K = !cdim ? 0 : K;
1658: (*T)->Nr = nrepl;
1659: (*T)->Np = npoints;
1660: (*T)->Nb = pdim;
1661: (*T)->Nc = Nc;
1662: (*T)->cdim = cdim;
1663: PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1664: for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1665: if (K >= 0) {
1666: for (p = 0; p < nrepl * npoints; ++p)
1667: for (d = 0; d < pdim; ++d)
1668: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1669: }
1670: if (K >= 1) {
1671: for (p = 0; p < nrepl * npoints; ++p)
1672: for (d = 0; d < pdim; ++d)
1673: for (c = 0; c < Nc; ++c)
1674: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1675: }
1676: if (K >= 2) {
1677: for (p = 0; p < nrepl * npoints; ++p)
1678: for (d = 0; d < pdim; ++d)
1679: for (c = 0; c < Nc; ++c)
1680: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1681: }
1682: PetscFunctionReturn(PETSC_SUCCESS);
1683: }
1685: /*@
1686: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1688: Input Parameters:
1689: + fvm - The `PetscFV` object
1690: . numFaces - The number of cell faces which are not constrained
1691: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1693: Output Parameter:
1694: . grad - the gradient
1696: Level: advanced
1698: .seealso: `PetscFV`, `PetscFVCreate()`
1699: @*/
1700: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1701: {
1702: PetscFunctionBegin;
1704: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1705: PetscFunctionReturn(PETSC_SUCCESS);
1706: }
1708: /*@C
1709: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1711: Not Collective
1713: Input Parameters:
1714: + fvm - The `PetscFV` object for the field being integrated
1715: . prob - The `PetscDS` specifying the discretizations and continuum functions
1716: . field - The field being integrated
1717: . Nf - The number of faces in the chunk
1718: . fgeom - The face geometry for each face in the chunk
1719: . neighborVol - The volume for each pair of cells in the chunk
1720: . uL - The state from the cell on the left
1721: - uR - The state from the cell on the right
1723: Output Parameters:
1724: + fluxL - the left fluxes for each face
1725: - fluxR - the right fluxes for each face
1727: Level: developer
1729: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1730: @*/
1731: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1732: {
1733: PetscFunctionBegin;
1735: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: /*@
1740: PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.
1742: Input Parameter:
1743: . fv - The initial `PetscFV`
1745: Output Parameter:
1746: . fvNew - A clone of the `PetscFV`
1748: Level: advanced
1750: Notes:
1751: This is typically used to change the number of components.
1753: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1754: @*/
1755: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1756: {
1757: PetscDualSpace Q;
1758: DM K;
1759: PetscQuadrature q;
1760: PetscInt Nc, cdim;
1762: PetscFunctionBegin;
1763: PetscCall(PetscFVGetDualSpace(fv, &Q));
1764: PetscCall(PetscFVGetQuadrature(fv, &q));
1765: PetscCall(PetscDualSpaceGetDM(Q, &K));
1767: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1768: PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1769: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1770: PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1771: PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1772: PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1773: PetscCall(PetscFVSetQuadrature(*fvNew, q));
1774: PetscFunctionReturn(PETSC_SUCCESS);
1775: }
1777: /*@
1778: PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1779: smaller copies.
1781: Input Parameter:
1782: . fv - The initial `PetscFV`
1784: Output Parameter:
1785: . fvRef - The refined `PetscFV`
1787: Level: advanced
1789: Notes:
1790: This is typically used to generate a preconditioner for a high order method from a lower order method on a
1791: refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1792: interpolation between regularly refined meshes.
1794: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1795: @*/
1796: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1797: {
1798: PetscDualSpace Q, Qref;
1799: DM K, Kref;
1800: PetscQuadrature q, qref;
1801: DMPolytopeType ct;
1802: DMPlexTransform tr;
1803: PetscReal *v0;
1804: PetscReal *jac, *invjac;
1805: PetscInt numComp, numSubelements, s;
1807: PetscFunctionBegin;
1808: PetscCall(PetscFVGetDualSpace(fv, &Q));
1809: PetscCall(PetscFVGetQuadrature(fv, &q));
1810: PetscCall(PetscDualSpaceGetDM(Q, &K));
1811: /* Create dual space */
1812: PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1813: PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1814: PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1815: PetscCall(DMDestroy(&Kref));
1816: PetscCall(PetscDualSpaceSetUp(Qref));
1817: /* Create volume */
1818: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1819: PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1820: PetscCall(PetscFVGetNumComponents(fv, &numComp));
1821: PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1822: PetscCall(PetscFVSetUp(*fvRef));
1823: /* Create quadrature */
1824: PetscCall(DMPlexGetCellType(K, 0, &ct));
1825: PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1826: PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1827: PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1828: PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1829: PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1830: for (s = 0; s < numSubelements; ++s) {
1831: PetscQuadrature qs;
1832: const PetscReal *points, *weights;
1833: PetscReal *p, *w;
1834: PetscInt dim, Nc, npoints, np;
1836: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1837: PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1838: np = npoints / numSubelements;
1839: PetscCall(PetscMalloc1(np * dim, &p));
1840: PetscCall(PetscMalloc1(np * Nc, &w));
1841: PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1842: PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1843: PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1844: PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1845: PetscCall(PetscQuadratureDestroy(&qs));
1846: }
1847: PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1848: PetscCall(DMPlexTransformDestroy(&tr));
1849: PetscCall(PetscQuadratureDestroy(&qref));
1850: PetscCall(PetscDualSpaceDestroy(&Qref));
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1855: {
1856: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1858: PetscFunctionBegin;
1859: PetscCall(PetscFree(b));
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1864: {
1865: PetscInt Nc, c;
1866: PetscViewerFormat format;
1868: PetscFunctionBegin;
1869: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1870: PetscCall(PetscViewerGetFormat(viewer, &format));
1871: PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1872: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1873: for (c = 0; c < Nc; c++) {
1874: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1875: }
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1880: {
1881: PetscBool iascii;
1883: PetscFunctionBegin;
1886: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1887: if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1888: PetscFunctionReturn(PETSC_SUCCESS);
1889: }
1891: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1892: {
1893: PetscFunctionBegin;
1894: PetscFunctionReturn(PETSC_SUCCESS);
1895: }
1897: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1898: {
1899: PetscInt dim;
1901: PetscFunctionBegin;
1902: PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1903: for (PetscInt f = 0; f < numFaces; ++f) {
1904: for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1905: }
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1909: /*
1910: neighborVol[f*2+0] contains the left geom
1911: neighborVol[f*2+1] contains the right geom
1912: */
1913: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1914: {
1915: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1916: void *rctx;
1917: PetscScalar *flux = fvm->fluxWork;
1918: const PetscScalar *constants;
1919: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1921: PetscFunctionBegin;
1922: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1923: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1924: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1925: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1926: PetscCall(PetscDSGetContext(prob, field, &rctx));
1927: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1928: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1929: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1930: for (f = 0; f < Nf; ++f) {
1931: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1932: for (d = 0; d < pdim; ++d) {
1933: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1934: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1935: }
1936: }
1937: PetscFunctionReturn(PETSC_SUCCESS);
1938: }
1940: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1941: {
1942: PetscFunctionBegin;
1943: fvm->ops->setfromoptions = NULL;
1944: fvm->ops->setup = PetscFVSetUp_Upwind;
1945: fvm->ops->view = PetscFVView_Upwind;
1946: fvm->ops->destroy = PetscFVDestroy_Upwind;
1947: fvm->ops->computegradient = PetscFVComputeGradient_Upwind;
1948: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: /*MC
1953: PETSCFVUPWIND = "upwind" - A `PetscFV` implementation
1955: Level: intermediate
1957: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1958: M*/
1960: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1961: {
1962: PetscFV_Upwind *b;
1964: PetscFunctionBegin;
1966: PetscCall(PetscNew(&b));
1967: fvm->data = b;
1969: PetscCall(PetscFVInitialize_Upwind(fvm));
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: #include <petscblaslapack.h>
1975: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1976: {
1977: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1979: PetscFunctionBegin;
1980: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1981: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1982: PetscCall(PetscFree(ls));
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1986: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1987: {
1988: PetscInt Nc, c;
1989: PetscViewerFormat format;
1991: PetscFunctionBegin;
1992: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1993: PetscCall(PetscViewerGetFormat(viewer, &format));
1994: PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1995: PetscCall(PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc));
1996: for (c = 0; c < Nc; c++) {
1997: if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1998: }
1999: PetscFunctionReturn(PETSC_SUCCESS);
2000: }
2002: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2003: {
2004: PetscBool iascii;
2006: PetscFunctionBegin;
2009: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2010: if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2011: PetscFunctionReturn(PETSC_SUCCESS);
2012: }
2014: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2015: {
2016: PetscFunctionBegin;
2017: PetscFunctionReturn(PETSC_SUCCESS);
2018: }
2020: /* Overwrites A. Can only handle full-rank problems with m>=n */
2021: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2022: {
2023: PetscBool debug = PETSC_FALSE;
2024: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2025: PetscScalar *R, *Q, *Aback, Alpha;
2027: PetscFunctionBegin;
2028: if (debug) {
2029: PetscCall(PetscMalloc1(m * n, &Aback));
2030: PetscCall(PetscArraycpy(Aback, A, m * n));
2031: }
2033: PetscCall(PetscBLASIntCast(m, &M));
2034: PetscCall(PetscBLASIntCast(n, &N));
2035: PetscCall(PetscBLASIntCast(mstride, &lda));
2036: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2037: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2038: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2039: PetscCall(PetscFPTrapPop());
2040: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2041: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2043: /* Extract an explicit representation of Q */
2044: Q = Ainv;
2045: PetscCall(PetscArraycpy(Q, A, mstride * n));
2046: K = N; /* full rank */
2047: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2048: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");
2050: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2051: Alpha = 1.0;
2052: ldb = lda;
2053: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2054: /* Ainv is Q, overwritten with inverse */
2056: if (debug) { /* Check that pseudo-inverse worked */
2057: PetscScalar Beta = 0.0;
2058: PetscBLASInt ldc;
2059: K = N;
2060: ldc = N;
2061: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2062: PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2063: PetscCall(PetscFree(Aback));
2064: }
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: /* Overwrites A. Can handle degenerate problems and m<n. */
2069: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2070: {
2071: PetscScalar *Brhs;
2072: PetscScalar *tmpwork;
2073: PetscReal rcond;
2074: #if defined(PETSC_USE_COMPLEX)
2075: PetscInt rworkSize;
2076: PetscReal *rwork;
2077: #endif
2078: PetscInt i, j, maxmn;
2079: PetscBLASInt M, N, lda, ldb, ldwork;
2080: PetscBLASInt nrhs, irank, info;
2082: PetscFunctionBegin;
2083: /* initialize to identity */
2084: tmpwork = work;
2085: Brhs = Ainv;
2086: maxmn = PetscMax(m, n);
2087: for (j = 0; j < maxmn; j++) {
2088: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2089: }
2091: PetscCall(PetscBLASIntCast(m, &M));
2092: PetscCall(PetscBLASIntCast(n, &N));
2093: PetscCall(PetscBLASIntCast(mstride, &lda));
2094: PetscCall(PetscBLASIntCast(maxmn, &ldb));
2095: PetscCall(PetscBLASIntCast(worksize, &ldwork));
2096: rcond = -1;
2097: nrhs = M;
2098: #if defined(PETSC_USE_COMPLEX)
2099: rworkSize = 5 * PetscMin(M, N);
2100: PetscCall(PetscMalloc1(rworkSize, &rwork));
2101: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2102: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2103: PetscCall(PetscFPTrapPop());
2104: PetscCall(PetscFree(rwork));
2105: #else
2106: nrhs = M;
2107: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2108: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2109: PetscCall(PetscFPTrapPop());
2110: #endif
2111: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2112: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2113: PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2114: PetscFunctionReturn(PETSC_SUCCESS);
2115: }
2117: #if 0
2118: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2119: {
2120: PetscReal grad[2] = {0, 0};
2121: const PetscInt *faces;
2122: PetscInt numFaces, f;
2124: PetscFunctionBegin;
2125: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2126: PetscCall(DMPlexGetCone(dm, cell, &faces));
2127: for (f = 0; f < numFaces; ++f) {
2128: const PetscInt *fcells;
2129: const CellGeom *cg1;
2130: const FaceGeom *fg;
2132: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2133: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2134: for (i = 0; i < 2; ++i) {
2135: PetscScalar du;
2137: if (fcells[i] == c) continue;
2138: PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2139: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2140: grad[0] += fg->grad[!i][0] * du;
2141: grad[1] += fg->grad[!i][1] * du;
2142: }
2143: }
2144: PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2147: #endif
2149: /*
2150: PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell
2152: Input Parameters:
2153: + fvm - The `PetscFV` object
2154: . numFaces - The number of cell faces which are not constrained
2155: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2157: Level: developer
2159: .seealso: `PetscFV`, `PetscFVCreate()`
2160: */
2161: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2162: {
2163: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2164: const PetscBool useSVD = PETSC_TRUE;
2165: const PetscInt maxFaces = ls->maxFaces;
2166: PetscInt dim, f, d;
2168: PetscFunctionBegin;
2169: if (numFaces > maxFaces) {
2170: PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2171: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2172: }
2173: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2174: for (f = 0; f < numFaces; ++f) {
2175: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2176: }
2177: /* Overwrites B with garbage, returns Binv in row-major format */
2178: if (useSVD) {
2179: PetscInt maxmn = PetscMax(numFaces, dim);
2180: PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2181: /* Binv shaped in column-major, coldim=maxmn.*/
2182: for (f = 0; f < numFaces; ++f) {
2183: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2184: }
2185: } else {
2186: PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2187: /* Binv shaped in row-major, rowdim=maxFaces.*/
2188: for (f = 0; f < numFaces; ++f) {
2189: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2190: }
2191: }
2192: PetscFunctionReturn(PETSC_SUCCESS);
2193: }
2195: /*
2196: neighborVol[f*2+0] contains the left geom
2197: neighborVol[f*2+1] contains the right geom
2198: */
2199: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2200: {
2201: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2202: void *rctx;
2203: PetscScalar *flux = fvm->fluxWork;
2204: const PetscScalar *constants;
2205: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2207: PetscFunctionBegin;
2208: PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2209: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2210: PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2211: PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2212: PetscCall(PetscDSGetContext(prob, field, &rctx));
2213: PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2214: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2215: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2216: for (f = 0; f < Nf; ++f) {
2217: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2218: for (d = 0; d < pdim; ++d) {
2219: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2220: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2221: }
2222: }
2223: PetscFunctionReturn(PETSC_SUCCESS);
2224: }
2226: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2227: {
2228: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2229: PetscInt dim, m, n, nrhs, minmn, maxmn;
2231: PetscFunctionBegin;
2233: PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2234: PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2235: ls->maxFaces = maxFaces;
2236: m = ls->maxFaces;
2237: n = dim;
2238: nrhs = ls->maxFaces;
2239: minmn = PetscMin(m, n);
2240: maxmn = PetscMax(m, n);
2241: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2242: PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2247: {
2248: PetscFunctionBegin;
2249: fvm->ops->setfromoptions = NULL;
2250: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2251: fvm->ops->view = PetscFVView_LeastSquares;
2252: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2253: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2254: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2255: PetscFunctionReturn(PETSC_SUCCESS);
2256: }
2258: /*MC
2259: PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation
2261: Level: intermediate
2263: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2264: M*/
2266: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2267: {
2268: PetscFV_LeastSquares *ls;
2270: PetscFunctionBegin;
2272: PetscCall(PetscNew(&ls));
2273: fvm->data = ls;
2275: ls->maxFaces = -1;
2276: ls->workSize = -1;
2277: ls->B = NULL;
2278: ls->Binv = NULL;
2279: ls->tau = NULL;
2280: ls->work = NULL;
2282: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2283: PetscCall(PetscFVInitialize_LeastSquares(fvm));
2284: PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2285: PetscFunctionReturn(PETSC_SUCCESS);
2286: }
2288: /*@
2289: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2291: Not Collective
2293: Input Parameters:
2294: + fvm - The `PetscFV` object
2295: - maxFaces - The maximum number of cell faces
2297: Level: intermediate
2299: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2300: @*/
2301: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2302: {
2303: PetscFunctionBegin;
2305: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2306: PetscFunctionReturn(PETSC_SUCCESS);
2307: }