Actual source code: dtfv.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/petscfvimpl.h> /*I "petscfv.h" I*/
2: #include <petsc/private/dmpleximpl.h> /* For CellRefiner */
3: #include <petscds.h>
5: PetscClassId PETSCLIMITER_CLASSID = 0;
7: PetscFunctionList PetscLimiterList = NULL;
8: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
10: PetscBool Limitercite = PETSC_FALSE;
11: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
12: " title = {Analysis of slope limiters on irregular grids},\n"
13: " journal = {AIAA paper},\n"
14: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
15: " volume = {490},\n"
16: " year = {2005}\n}\n";
20: /*@C
21: PetscLimiterRegister - Adds a new PetscLimiter implementation
23: Not Collective
25: Input Parameters:
26: + name - The name of a new user-defined creation routine
27: - create_func - The creation routine itself
29: Notes:
30: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
32: Sample usage:
33: .vb
34: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
35: .ve
37: Then, your PetscLimiter type can be chosen with the procedural interface via
38: .vb
39: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
40: PetscLimiterSetType(PetscLimiter, "my_lim");
41: .ve
42: or at runtime via the option
43: .vb
44: -petsclimiter_type my_lim
45: .ve
47: Level: advanced
49: .keywords: PetscLimiter, register
50: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
52: @*/
53: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
54: {
58: PetscFunctionListAdd(&PetscLimiterList, sname, function);
59: return(0);
60: }
64: /*@C
65: PetscLimiterSetType - Builds a particular PetscLimiter
67: Collective on PetscLimiter
69: Input Parameters:
70: + lim - The PetscLimiter object
71: - name - The kind of limiter
73: Options Database Key:
74: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
76: Level: intermediate
78: .keywords: PetscLimiter, set, type
79: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
80: @*/
81: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
82: {
83: PetscErrorCode (*r)(PetscLimiter);
84: PetscBool match;
89: PetscObjectTypeCompare((PetscObject) lim, name, &match);
90: if (match) return(0);
92: PetscLimiterRegisterAll();
93: PetscFunctionListFind(PetscLimiterList, name, &r);
94: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
96: if (lim->ops->destroy) {
97: (*lim->ops->destroy)(lim);
98: lim->ops->destroy = NULL;
99: }
100: (*r)(lim);
101: PetscObjectChangeTypeName((PetscObject) lim, name);
102: return(0);
103: }
107: /*@C
108: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
110: Not Collective
112: Input Parameter:
113: . lim - The PetscLimiter
115: Output Parameter:
116: . name - The PetscLimiter type name
118: Level: intermediate
120: .keywords: PetscLimiter, get, type, name
121: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
122: @*/
123: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
124: {
130: PetscLimiterRegisterAll();
131: *name = ((PetscObject) lim)->type_name;
132: return(0);
133: }
137: /*@C
138: PetscLimiterView - Views a PetscLimiter
140: Collective on PetscLimiter
142: Input Parameter:
143: + lim - the PetscLimiter object to view
144: - v - the viewer
146: Level: developer
148: .seealso: PetscLimiterDestroy()
149: @*/
150: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
151: {
156: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
157: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
158: return(0);
159: }
163: /*@
164: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
166: Collective on PetscLimiter
168: Input Parameter:
169: . lim - the PetscLimiter object to set options for
171: Level: developer
173: .seealso: PetscLimiterView()
174: @*/
175: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
176: {
177: const char *defaultType;
178: char name[256];
179: PetscBool flg;
184: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
185: else defaultType = ((PetscObject) lim)->type_name;
186: PetscLimiterRegisterAll();
188: PetscObjectOptionsBegin((PetscObject) lim);
189: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
190: if (flg) {
191: PetscLimiterSetType(lim, name);
192: } else if (!((PetscObject) lim)->type_name) {
193: PetscLimiterSetType(lim, defaultType);
194: }
195: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
196: /* process any options handlers added with PetscObjectAddOptionsHandler() */
197: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
198: PetscOptionsEnd();
199: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
200: return(0);
201: }
205: /*@C
206: PetscLimiterSetUp - Construct data structures for the PetscLimiter
208: Collective on PetscLimiter
210: Input Parameter:
211: . lim - the PetscLimiter object to setup
213: Level: developer
215: .seealso: PetscLimiterView(), PetscLimiterDestroy()
216: @*/
217: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
218: {
223: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
224: return(0);
225: }
229: /*@
230: PetscLimiterDestroy - Destroys a PetscLimiter object
232: Collective on PetscLimiter
234: Input Parameter:
235: . lim - the PetscLimiter object to destroy
237: Level: developer
239: .seealso: PetscLimiterView()
240: @*/
241: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
242: {
246: if (!*lim) return(0);
249: if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
250: ((PetscObject) (*lim))->refct = 0;
252: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
253: PetscHeaderDestroy(lim);
254: return(0);
255: }
259: /*@
260: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
262: Collective on MPI_Comm
264: Input Parameter:
265: . comm - The communicator for the PetscLimiter object
267: Output Parameter:
268: . lim - The PetscLimiter object
270: Level: beginner
272: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
273: @*/
274: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
275: {
276: PetscLimiter l;
281: PetscCitationsRegister(LimiterCitation,&Limitercite);
282: *lim = NULL;
283: PetscFVInitializePackage();
285: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
287: *lim = l;
288: return(0);
289: }
293: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
294: *
295: * The classical flux-limited formulation is psi(r) where
296: *
297: * r = (u[0] - u[-1]) / (u[1] - u[0])
298: *
299: * The second order TVD region is bounded by
300: *
301: * psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
302: *
303: * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
304: * phi(r)(r+1)/2 in which the reconstructed interface values are
305: *
306: * u(v) = u[0] + phi(r) (grad u)[0] v
307: *
308: * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
309: *
310: * 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))
311: *
312: * For a nicer symmetric formulation, rewrite in terms of
313: *
314: * f = (u[0] - u[-1]) / (u[1] - u[-1])
315: *
316: * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
317: *
318: * phi(r) = phi(1/r)
319: *
320: * becomes
321: *
322: * w(f) = w(1-f).
323: *
324: * The limiters below implement this final form w(f). The reference methods are
325: *
326: * w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
327: * */
328: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
329: {
335: (*lim->ops->limit)(lim, flim, phi);
336: return(0);
337: }
341: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
342: {
343: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
344: PetscErrorCode ierr;
347: PetscFree(l);
348: return(0);
349: }
353: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
354: {
355: PetscViewerFormat format;
356: PetscErrorCode ierr;
359: PetscViewerGetFormat(viewer, &format);
360: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
361: return(0);
362: }
366: PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
367: {
368: PetscBool iascii;
374: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
375: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
376: return(0);
377: }
381: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
382: {
384: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
385: return(0);
386: }
390: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
391: {
393: lim->ops->view = PetscLimiterView_Sin;
394: lim->ops->destroy = PetscLimiterDestroy_Sin;
395: lim->ops->limit = PetscLimiterLimit_Sin;
396: return(0);
397: }
399: /*MC
400: PETSCLIMITERSIN = "sin" - A PetscLimiter object
402: Level: intermediate
404: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
405: M*/
409: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
410: {
411: PetscLimiter_Sin *l;
412: PetscErrorCode ierr;
416: PetscNewLog(lim, &l);
417: lim->data = l;
419: PetscLimiterInitialize_Sin(lim);
420: return(0);
421: }
425: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
426: {
427: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
428: PetscErrorCode ierr;
431: PetscFree(l);
432: return(0);
433: }
437: PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
438: {
439: PetscViewerFormat format;
440: PetscErrorCode ierr;
443: PetscViewerGetFormat(viewer, &format);
444: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
445: return(0);
446: }
450: PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
451: {
452: PetscBool iascii;
458: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
459: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
460: return(0);
461: }
465: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
466: {
468: *phi = 0.0;
469: return(0);
470: }
474: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
475: {
477: lim->ops->view = PetscLimiterView_Zero;
478: lim->ops->destroy = PetscLimiterDestroy_Zero;
479: lim->ops->limit = PetscLimiterLimit_Zero;
480: return(0);
481: }
483: /*MC
484: PETSCLIMITERZERO = "zero" - A PetscLimiter object
486: Level: intermediate
488: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
489: M*/
493: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
494: {
495: PetscLimiter_Zero *l;
496: PetscErrorCode ierr;
500: PetscNewLog(lim, &l);
501: lim->data = l;
503: PetscLimiterInitialize_Zero(lim);
504: return(0);
505: }
509: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
510: {
511: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
512: PetscErrorCode ierr;
515: PetscFree(l);
516: return(0);
517: }
521: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
522: {
523: PetscViewerFormat format;
524: PetscErrorCode ierr;
527: PetscViewerGetFormat(viewer, &format);
528: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
529: return(0);
530: }
534: PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
535: {
536: PetscBool iascii;
542: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
543: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
544: return(0);
545: }
549: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
550: {
552: *phi = 1.0;
553: return(0);
554: }
558: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
559: {
561: lim->ops->view = PetscLimiterView_None;
562: lim->ops->destroy = PetscLimiterDestroy_None;
563: lim->ops->limit = PetscLimiterLimit_None;
564: return(0);
565: }
567: /*MC
568: PETSCLIMITERNONE = "none" - A PetscLimiter object
570: Level: intermediate
572: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
573: M*/
577: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
578: {
579: PetscLimiter_None *l;
580: PetscErrorCode ierr;
584: PetscNewLog(lim, &l);
585: lim->data = l;
587: PetscLimiterInitialize_None(lim);
588: return(0);
589: }
593: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
594: {
595: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
596: PetscErrorCode ierr;
599: PetscFree(l);
600: return(0);
601: }
605: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
606: {
607: PetscViewerFormat format;
608: PetscErrorCode ierr;
611: PetscViewerGetFormat(viewer, &format);
612: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
613: return(0);
614: }
618: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
619: {
620: PetscBool iascii;
626: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
627: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
628: return(0);
629: }
633: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
634: {
636: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
637: return(0);
638: }
642: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
643: {
645: lim->ops->view = PetscLimiterView_Minmod;
646: lim->ops->destroy = PetscLimiterDestroy_Minmod;
647: lim->ops->limit = PetscLimiterLimit_Minmod;
648: return(0);
649: }
651: /*MC
652: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
654: Level: intermediate
656: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
657: M*/
661: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
662: {
663: PetscLimiter_Minmod *l;
664: PetscErrorCode ierr;
668: PetscNewLog(lim, &l);
669: lim->data = l;
671: PetscLimiterInitialize_Minmod(lim);
672: return(0);
673: }
677: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
678: {
679: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
680: PetscErrorCode ierr;
683: PetscFree(l);
684: return(0);
685: }
689: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
690: {
691: PetscViewerFormat format;
692: PetscErrorCode ierr;
695: PetscViewerGetFormat(viewer, &format);
696: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
697: return(0);
698: }
702: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
703: {
704: PetscBool iascii;
710: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
711: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
712: return(0);
713: }
717: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
718: {
720: *phi = PetscMax(0, 4*f*(1-f));
721: return(0);
722: }
726: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
727: {
729: lim->ops->view = PetscLimiterView_VanLeer;
730: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
731: lim->ops->limit = PetscLimiterLimit_VanLeer;
732: return(0);
733: }
735: /*MC
736: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
738: Level: intermediate
740: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
741: M*/
745: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
746: {
747: PetscLimiter_VanLeer *l;
748: PetscErrorCode ierr;
752: PetscNewLog(lim, &l);
753: lim->data = l;
755: PetscLimiterInitialize_VanLeer(lim);
756: return(0);
757: }
761: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
762: {
763: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
764: PetscErrorCode ierr;
767: PetscFree(l);
768: return(0);
769: }
773: PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
774: {
775: PetscViewerFormat format;
776: PetscErrorCode ierr;
779: PetscViewerGetFormat(viewer, &format);
780: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
781: return(0);
782: }
786: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
787: {
788: PetscBool iascii;
794: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
795: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
796: return(0);
797: }
801: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
802: {
804: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
805: return(0);
806: }
810: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
811: {
813: lim->ops->view = PetscLimiterView_VanAlbada;
814: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
815: lim->ops->limit = PetscLimiterLimit_VanAlbada;
816: return(0);
817: }
819: /*MC
820: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
822: Level: intermediate
824: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
825: M*/
829: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
830: {
831: PetscLimiter_VanAlbada *l;
832: PetscErrorCode ierr;
836: PetscNewLog(lim, &l);
837: lim->data = l;
839: PetscLimiterInitialize_VanAlbada(lim);
840: return(0);
841: }
845: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
846: {
847: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
848: PetscErrorCode ierr;
851: PetscFree(l);
852: return(0);
853: }
857: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
858: {
859: PetscViewerFormat format;
860: PetscErrorCode ierr;
863: PetscViewerGetFormat(viewer, &format);
864: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
865: return(0);
866: }
870: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
871: {
872: PetscBool iascii;
878: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
879: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
880: return(0);
881: }
885: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
886: {
888: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
889: return(0);
890: }
894: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
895: {
897: lim->ops->view = PetscLimiterView_Superbee;
898: lim->ops->destroy = PetscLimiterDestroy_Superbee;
899: lim->ops->limit = PetscLimiterLimit_Superbee;
900: return(0);
901: }
903: /*MC
904: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
906: Level: intermediate
908: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
909: M*/
913: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
914: {
915: PetscLimiter_Superbee *l;
916: PetscErrorCode ierr;
920: PetscNewLog(lim, &l);
921: lim->data = l;
923: PetscLimiterInitialize_Superbee(lim);
924: return(0);
925: }
929: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
930: {
931: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
932: PetscErrorCode ierr;
935: PetscFree(l);
936: return(0);
937: }
941: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
942: {
943: PetscViewerFormat format;
944: PetscErrorCode ierr;
947: PetscViewerGetFormat(viewer, &format);
948: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
949: return(0);
950: }
954: PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
955: {
956: PetscBool iascii;
962: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
963: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
964: return(0);
965: }
969: /* aka Barth-Jespersen */
970: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
971: {
973: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
974: return(0);
975: }
979: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
980: {
982: lim->ops->view = PetscLimiterView_MC;
983: lim->ops->destroy = PetscLimiterDestroy_MC;
984: lim->ops->limit = PetscLimiterLimit_MC;
985: return(0);
986: }
988: /*MC
989: PETSCLIMITERMC = "mc" - A PetscLimiter object
991: Level: intermediate
993: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
994: M*/
998: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
999: {
1000: PetscLimiter_MC *l;
1001: PetscErrorCode ierr;
1005: PetscNewLog(lim, &l);
1006: lim->data = l;
1008: PetscLimiterInitialize_MC(lim);
1009: return(0);
1010: }
1012: PetscClassId PETSCFV_CLASSID = 0;
1014: PetscFunctionList PetscFVList = NULL;
1015: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
1019: /*@C
1020: PetscFVRegister - Adds a new PetscFV implementation
1022: Not Collective
1024: Input Parameters:
1025: + name - The name of a new user-defined creation routine
1026: - create_func - The creation routine itself
1028: Notes:
1029: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
1031: Sample usage:
1032: .vb
1033: PetscFVRegister("my_fv", MyPetscFVCreate);
1034: .ve
1036: Then, your PetscFV type can be chosen with the procedural interface via
1037: .vb
1038: PetscFVCreate(MPI_Comm, PetscFV *);
1039: PetscFVSetType(PetscFV, "my_fv");
1040: .ve
1041: or at runtime via the option
1042: .vb
1043: -petscfv_type my_fv
1044: .ve
1046: Level: advanced
1048: .keywords: PetscFV, register
1049: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
1051: @*/
1052: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
1053: {
1057: PetscFunctionListAdd(&PetscFVList, sname, function);
1058: return(0);
1059: }
1063: /*@C
1064: PetscFVSetType - Builds a particular PetscFV
1066: Collective on PetscFV
1068: Input Parameters:
1069: + fvm - The PetscFV object
1070: - name - The kind of FVM space
1072: Options Database Key:
1073: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
1075: Level: intermediate
1077: .keywords: PetscFV, set, type
1078: .seealso: PetscFVGetType(), PetscFVCreate()
1079: @*/
1080: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
1081: {
1082: PetscErrorCode (*r)(PetscFV);
1083: PetscBool match;
1088: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1089: if (match) return(0);
1091: PetscFVRegisterAll();
1092: PetscFunctionListFind(PetscFVList, name, &r);
1093: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1095: if (fvm->ops->destroy) {
1096: (*fvm->ops->destroy)(fvm);
1097: fvm->ops->destroy = NULL;
1098: }
1099: (*r)(fvm);
1100: PetscObjectChangeTypeName((PetscObject) fvm, name);
1101: return(0);
1102: }
1106: /*@C
1107: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1109: Not Collective
1111: Input Parameter:
1112: . fvm - The PetscFV
1114: Output Parameter:
1115: . name - The PetscFV type name
1117: Level: intermediate
1119: .keywords: PetscFV, get, type, name
1120: .seealso: PetscFVSetType(), PetscFVCreate()
1121: @*/
1122: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1123: {
1129: PetscFVRegisterAll();
1130: *name = ((PetscObject) fvm)->type_name;
1131: return(0);
1132: }
1136: /*@C
1137: PetscFVView - Views a PetscFV
1139: Collective on PetscFV
1141: Input Parameter:
1142: + fvm - the PetscFV object to view
1143: - v - the viewer
1145: Level: developer
1147: .seealso: PetscFVDestroy()
1148: @*/
1149: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1150: {
1155: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1156: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1157: return(0);
1158: }
1162: /*@
1163: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1165: Collective on PetscFV
1167: Input Parameter:
1168: . fvm - the PetscFV object to set options for
1170: Level: developer
1172: .seealso: PetscFVView()
1173: @*/
1174: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1175: {
1176: const char *defaultType;
1177: char name[256];
1178: PetscBool flg;
1183: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1184: else defaultType = ((PetscObject) fvm)->type_name;
1185: PetscFVRegisterAll();
1187: PetscObjectOptionsBegin((PetscObject) fvm);
1188: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1189: if (flg) {
1190: PetscFVSetType(fvm, name);
1191: } else if (!((PetscObject) fvm)->type_name) {
1192: PetscFVSetType(fvm, defaultType);
1193: }
1194: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1195: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1196: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1197: PetscLimiterSetFromOptions(fvm->limiter);
1198: PetscOptionsEnd();
1199: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1200: return(0);
1201: }
1205: /*@
1206: PetscFVSetUp - Construct data structures for the PetscFV
1208: Collective on PetscFV
1210: Input Parameter:
1211: . fvm - the PetscFV object to setup
1213: Level: developer
1215: .seealso: PetscFVView(), PetscFVDestroy()
1216: @*/
1217: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1218: {
1223: PetscLimiterSetUp(fvm->limiter);
1224: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1225: return(0);
1226: }
1230: /*@
1231: PetscFVDestroy - Destroys a PetscFV object
1233: Collective on PetscFV
1235: Input Parameter:
1236: . fvm - the PetscFV object to destroy
1238: Level: developer
1240: .seealso: PetscFVView()
1241: @*/
1242: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1243: {
1244: PetscInt i;
1248: if (!*fvm) return(0);
1251: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1252: ((PetscObject) (*fvm))->refct = 0;
1254: for (i = 0; i < (*fvm)->numComponents; i++) {
1255: PetscFree((*fvm)->componentNames[i]);
1256: }
1257: PetscFree((*fvm)->componentNames);
1258: PetscLimiterDestroy(&(*fvm)->limiter);
1259: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1260: PetscFree((*fvm)->fluxWork);
1261: PetscQuadratureDestroy(&(*fvm)->quadrature);
1262: PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);
1264: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1265: PetscHeaderDestroy(fvm);
1266: return(0);
1267: }
1271: /*@
1272: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1274: Collective on MPI_Comm
1276: Input Parameter:
1277: . comm - The communicator for the PetscFV object
1279: Output Parameter:
1280: . fvm - The PetscFV object
1282: Level: beginner
1284: .seealso: PetscFVSetType(), PETSCFVUPWIND
1285: @*/
1286: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1287: {
1288: PetscFV f;
1293: *fvm = NULL;
1294: PetscFVInitializePackage();
1296: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1297: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1299: PetscLimiterCreate(comm, &f->limiter);
1300: f->numComponents = 1;
1301: f->dim = 0;
1302: f->computeGradients = PETSC_FALSE;
1303: f->fluxWork = NULL;
1304: PetscCalloc1(f->numComponents,&f->componentNames);
1306: *fvm = f;
1307: return(0);
1308: }
1312: /*@
1313: PetscFVSetLimiter - Set the limiter object
1315: Logically collective on PetscFV
1317: Input Parameters:
1318: + fvm - the PetscFV object
1319: - lim - The PetscLimiter
1321: Level: developer
1323: .seealso: PetscFVGetLimiter()
1324: @*/
1325: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1326: {
1332: PetscLimiterDestroy(&fvm->limiter);
1333: PetscObjectReference((PetscObject) lim);
1334: fvm->limiter = lim;
1335: return(0);
1336: }
1340: /*@
1341: PetscFVGetLimiter - Get the limiter object
1343: Not collective
1345: Input Parameter:
1346: . fvm - the PetscFV object
1348: Output Parameter:
1349: . lim - The PetscLimiter
1351: Level: developer
1353: .seealso: PetscFVSetLimiter()
1354: @*/
1355: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1356: {
1360: *lim = fvm->limiter;
1361: return(0);
1362: }
1366: /*@
1367: PetscFVSetNumComponents - Set the number of field components
1369: Logically collective on PetscFV
1371: Input Parameters:
1372: + fvm - the PetscFV object
1373: - comp - The number of components
1375: Level: developer
1377: .seealso: PetscFVGetNumComponents()
1378: @*/
1379: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1380: {
1385: if (fvm->numComponents != comp) {
1386: PetscInt i;
1388: for (i = 0; i < fvm->numComponents; i++) {
1389: PetscFree(fvm->componentNames[i]);
1390: }
1391: PetscFree(fvm->componentNames);
1392: PetscCalloc1(comp,&fvm->componentNames);
1393: }
1394: fvm->numComponents = comp;
1395: PetscFree(fvm->fluxWork);
1396: PetscMalloc1(comp, &fvm->fluxWork);
1397: return(0);
1398: }
1402: /*@
1403: PetscFVGetNumComponents - Get the number of field components
1405: Not collective
1407: Input Parameter:
1408: . fvm - the PetscFV object
1410: Output Parameter:
1411: , comp - The number of components
1413: Level: developer
1415: .seealso: PetscFVSetNumComponents()
1416: @*/
1417: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1418: {
1422: *comp = fvm->numComponents;
1423: return(0);
1424: }
1428: /*@C
1429: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1431: Logically collective on PetscFV
1432: Input Parameters:
1433: + fvm - the PetscFV object
1434: . comp - the component number
1435: - name - the component name
1437: Level: developer
1439: .seealso: PetscFVGetComponentName()
1440: @*/
1441: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1442: {
1446: PetscFree(fvm->componentNames[comp]);
1447: PetscStrallocpy(name,&fvm->componentNames[comp]);
1448: return(0);
1449: }
1453: /*@C
1454: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1456: Logically collective on PetscFV
1457: Input Parameters:
1458: + fvm - the PetscFV object
1459: - comp - the component number
1461: Output Parameter:
1462: . name - the component name
1464: Level: developer
1466: .seealso: PetscFVSetComponentName()
1467: @*/
1468: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1469: {
1471: *name = fvm->componentNames[comp];
1472: return(0);
1473: }
1477: /*@
1478: PetscFVSetSpatialDimension - Set the spatial dimension
1480: Logically collective on PetscFV
1482: Input Parameters:
1483: + fvm - the PetscFV object
1484: - dim - The spatial dimension
1486: Level: developer
1488: .seealso: PetscFVGetSpatialDimension()
1489: @*/
1490: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1491: {
1494: fvm->dim = dim;
1495: return(0);
1496: }
1500: /*@
1501: PetscFVGetSpatialDimension - Get the spatial dimension
1503: Logically collective on PetscFV
1505: Input Parameter:
1506: . fvm - the PetscFV object
1508: Output Parameter:
1509: . dim - The spatial dimension
1511: Level: developer
1513: .seealso: PetscFVSetSpatialDimension()
1514: @*/
1515: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1516: {
1520: *dim = fvm->dim;
1521: return(0);
1522: }
1526: /*@
1527: PetscFVSetComputeGradients - Toggle computation of cell gradients
1529: Logically collective on PetscFV
1531: Input Parameters:
1532: + fvm - the PetscFV object
1533: - computeGradients - Flag to compute cell gradients
1535: Level: developer
1537: .seealso: PetscFVGetComputeGradients()
1538: @*/
1539: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1540: {
1543: fvm->computeGradients = computeGradients;
1544: return(0);
1545: }
1549: /*@
1550: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1552: Not collective
1554: Input Parameter:
1555: . fvm - the PetscFV object
1557: Output Parameter:
1558: . computeGradients - Flag to compute cell gradients
1560: Level: developer
1562: .seealso: PetscFVSetComputeGradients()
1563: @*/
1564: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1565: {
1569: *computeGradients = fvm->computeGradients;
1570: return(0);
1571: }
1575: /*@
1576: PetscFVSetQuadrature - Set the quadrature object
1578: Logically collective on PetscFV
1580: Input Parameters:
1581: + fvm - the PetscFV object
1582: - q - The PetscQuadrature
1584: Level: developer
1586: .seealso: PetscFVGetQuadrature()
1587: @*/
1588: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1589: {
1594: PetscQuadratureDestroy(&fvm->quadrature);
1595: PetscObjectReference((PetscObject) q);
1596: fvm->quadrature = q;
1597: return(0);
1598: }
1602: /*@
1603: PetscFVGetQuadrature - Get the quadrature object
1605: Not collective
1607: Input Parameter:
1608: . fvm - the PetscFV object
1610: Output Parameter:
1611: . lim - The PetscQuadrature
1613: Level: developer
1615: .seealso: PetscFVSetQuadrature()
1616: @*/
1617: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1618: {
1622: if (!fvm->quadrature) {
1623: /* Create default 1-point quadrature */
1624: PetscReal *points, *weights;
1627: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1628: PetscCalloc1(fvm->dim, &points);
1629: PetscMalloc1(1, &weights);
1630: weights[0] = 1.0;
1631: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, points, weights);
1632: }
1633: *q = fvm->quadrature;
1634: return(0);
1635: }
1639: /*@
1640: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1642: Not collective
1644: Input Parameter:
1645: . fvm - The PetscFV object
1647: Output Parameter:
1648: . sp - The PetscDualSpace object
1650: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1652: Level: developer
1654: .seealso: PetscFVCreate()
1655: @*/
1656: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1657: {
1661: if (!fvm->dualSpace) {
1662: DM K;
1663: PetscQuadrature q;
1664: PetscInt dim;
1665: PetscErrorCode ierr;
1667: PetscFVGetSpatialDimension(fvm, &dim);
1668: PetscFVGetQuadrature(fvm, &q);
1669: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1670: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1671: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1672: PetscDualSpaceSetDM(fvm->dualSpace, K);
1673: DMDestroy(&K);
1674: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, 1);
1675: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, 0, q);
1676: }
1677: *sp = fvm->dualSpace;
1678: return(0);
1679: }
1683: /*@
1684: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1686: Not collective
1688: Input Parameters:
1689: + fvm - The PetscFV object
1690: - sp - The PetscDualSpace object
1692: Level: intermediate
1694: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1696: .seealso: PetscFVCreate()
1697: @*/
1698: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1699: {
1705: PetscDualSpaceDestroy(&fvm->dualSpace);
1706: fvm->dualSpace = sp;
1707: PetscObjectReference((PetscObject) fvm->dualSpace);
1708: return(0);
1709: }
1713: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1714: {
1715: PetscInt npoints;
1716: const PetscReal *points;
1717: PetscErrorCode ierr;
1724: PetscQuadratureGetData(fvm->quadrature, NULL, &npoints, &points, NULL);
1725: if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1726: if (B) *B = fvm->B;
1727: if (D) *D = fvm->D;
1728: if (H) *H = fvm->H;
1729: return(0);
1730: }
1734: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1735: {
1736: PetscInt pdim = 1; /* Dimension of approximation space P */
1737: PetscInt dim; /* Spatial dimension */
1738: PetscInt comp; /* Field components */
1739: PetscInt p, d, c, e;
1740: PetscErrorCode ierr;
1748: PetscFVGetSpatialDimension(fvm, &dim);
1749: PetscFVGetNumComponents(fvm, &comp);
1750: if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1751: if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1752: if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1753: if (B) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) (*B)[(p*pdim + d)*comp + c] = 1.0;}
1754: if (D) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim; ++e) (*D)[((p*pdim + d)*comp + c)*dim + e] = 0.0;}
1755: if (H) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim*dim; ++e) (*H)[((p*pdim + d)*comp + c)*dim*dim + e] = 0.0;}
1756: return(0);
1757: }
1761: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1762: {
1767: if (B && *B) {PetscFree(*B);}
1768: if (D && *D) {PetscFree(*D);}
1769: if (H && *H) {PetscFree(*H);}
1770: return(0);
1771: }
1775: /*@C
1776: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1778: Input Parameters:
1779: + fvm - The PetscFV object
1780: . numFaces - The number of cell faces which are not constrained
1781: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1783: Level: developer
1785: .seealso: PetscFVCreate()
1786: @*/
1787: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1788: {
1793: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1794: return(0);
1795: }
1799: /*C
1800: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1802: Not collective
1804: Input Parameters:
1805: + fvm - The PetscFV object for the field being integrated
1806: . prob - The PetscDS specifing the discretizations and continuum functions
1807: . field - The field being integrated
1808: . Nf - The number of faces in the chunk
1809: . fgeom - The face geometry for each face in the chunk
1810: . neighborVol - The volume for each pair of cells in the chunk
1811: . uL - The state from the cell on the left
1812: - uR - The state from the cell on the right
1814: Output Parameter
1815: + fluxL - the left fluxes for each face
1816: - fluxR - the right fluxes for each face
1817: */
1818: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1819: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1820: {
1825: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1826: return(0);
1827: }
1831: /*@
1832: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1833: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1834: sparsity). It is also used to create an interpolation between regularly refined meshes.
1836: Input Parameter:
1837: . fv - The initial PetscFV
1839: Output Parameter:
1840: . fvRef - The refined PetscFV
1842: Level: developer
1844: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1845: @*/
1846: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1847: {
1848: PetscDualSpace Q, Qref;
1849: DM K, Kref;
1850: PetscQuadrature q, qref;
1851: CellRefiner cellRefiner;
1852: PetscReal *v0;
1853: PetscReal *jac, *invjac;
1854: PetscInt numComp, numSubelements, s;
1855: PetscErrorCode ierr;
1858: PetscFVGetDualSpace(fv, &Q);
1859: PetscFVGetQuadrature(fv, &q);
1860: PetscDualSpaceGetDM(Q, &K);
1861: /* Create dual space */
1862: PetscDualSpaceDuplicate(Q, &Qref);
1863: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1864: PetscDualSpaceSetDM(Qref, Kref);
1865: DMDestroy(&Kref);
1866: PetscDualSpaceSetUp(Qref);
1867: /* Create volume */
1868: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1869: PetscFVSetDualSpace(*fvRef, Qref);
1870: PetscFVGetNumComponents(fv, &numComp);
1871: PetscFVSetNumComponents(*fvRef, numComp);
1872: PetscFVSetUp(*fvRef);
1873: /* Create quadrature */
1874: DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1875: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1876: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1877: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1878: for (s = 0; s < numSubelements; ++s) {
1879: PetscQuadrature qs;
1880: const PetscReal *points, *weights;
1881: PetscReal *p, *w;
1882: PetscInt dim, npoints, np;
1884: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1885: PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
1886: np = npoints/numSubelements;
1887: PetscMalloc1(np*dim,&p);
1888: PetscMalloc1(np,&w);
1889: PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1890: PetscMemcpy(w, &weights[s*np], np * sizeof(PetscReal));
1891: PetscQuadratureSetData(qs, dim, np, p, w);
1892: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1893: PetscQuadratureDestroy(&qs);
1894: }
1895: CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1896: PetscFVSetQuadrature(*fvRef, qref);
1897: PetscQuadratureDestroy(&qref);
1898: PetscDualSpaceDestroy(&Qref);
1899: return(0);
1900: }
1904: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1905: {
1906: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1910: PetscFree(b);
1911: return(0);
1912: }
1916: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1917: {
1918: PetscInt Nc, c;
1919: PetscViewerFormat format;
1920: PetscErrorCode ierr;
1923: PetscFVGetNumComponents(fv, &Nc);
1924: PetscViewerGetFormat(viewer, &format);
1925: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1926: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1927: for (c = 0; c < Nc; c++) {
1928: if (fv->componentNames[c]) {
1929: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1930: }
1931: }
1932: return(0);
1933: }
1937: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1938: {
1939: PetscBool iascii;
1945: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1946: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1947: return(0);
1948: }
1952: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1953: {
1955: return(0);
1956: }
1960: /*
1961: neighborVol[f*2+0] contains the left geom
1962: neighborVol[f*2+1] contains the right geom
1963: */
1964: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1965: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1966: {
1967: void (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
1968: void *rctx;
1969: PetscScalar *flux = fvm->fluxWork;
1970: PetscInt dim, pdim, totDim, Nc, off, f, d;
1974: PetscDSGetTotalComponents(prob, &Nc);
1975: PetscDSGetTotalDimension(prob, &totDim);
1976: PetscDSGetFieldOffset(prob, field, &off);
1977: PetscDSGetRiemannSolver(prob, field, &riemann);
1978: PetscDSGetContext(prob, field, &rctx);
1979: PetscFVGetSpatialDimension(fvm, &dim);
1980: PetscFVGetNumComponents(fvm, &pdim);
1981: for (f = 0; f < Nf; ++f) {
1982: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
1983: for (d = 0; d < pdim; ++d) {
1984: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1985: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1986: }
1987: }
1988: return(0);
1989: }
1993: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1994: {
1996: fvm->ops->setfromoptions = NULL;
1997: fvm->ops->setup = PetscFVSetUp_Upwind;
1998: fvm->ops->view = PetscFVView_Upwind;
1999: fvm->ops->destroy = PetscFVDestroy_Upwind;
2000: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
2001: return(0);
2002: }
2004: /*MC
2005: PETSCFVUPWIND = "upwind" - A PetscFV object
2007: Level: intermediate
2009: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2010: M*/
2014: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
2015: {
2016: PetscFV_Upwind *b;
2017: PetscErrorCode ierr;
2021: PetscNewLog(fvm,&b);
2022: fvm->data = b;
2024: PetscFVInitialize_Upwind(fvm);
2025: return(0);
2026: }
2028: #include <petscblaslapack.h>
2032: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
2033: {
2034: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2035: PetscErrorCode ierr;
2038: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
2039: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2040: PetscFree(ls);
2041: return(0);
2042: }
2046: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
2047: {
2048: PetscInt Nc, c;
2049: PetscViewerFormat format;
2050: PetscErrorCode ierr;
2053: PetscFVGetNumComponents(fv, &Nc);
2054: PetscViewerGetFormat(viewer, &format);
2055: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
2056: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
2057: for (c = 0; c < Nc; c++) {
2058: if (fv->componentNames[c]) {
2059: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
2060: }
2061: }
2062: return(0);
2063: }
2067: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2068: {
2069: PetscBool iascii;
2075: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2076: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2077: return(0);
2078: }
2082: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2083: {
2085: return(0);
2086: }
2090: /* Overwrites A. Can only handle full-rank problems with m>=n */
2091: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2092: {
2093: PetscBool debug = PETSC_FALSE;
2095: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2096: PetscScalar *R,*Q,*Aback,Alpha;
2099: if (debug) {
2100: PetscMalloc1(m*n,&Aback);
2101: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2102: }
2104: PetscBLASIntCast(m,&M);
2105: PetscBLASIntCast(n,&N);
2106: PetscBLASIntCast(mstride,&lda);
2107: PetscBLASIntCast(worksize,&ldwork);
2108: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2109: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2110: PetscFPTrapPop();
2111: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2112: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2114: /* Extract an explicit representation of Q */
2115: Q = Ainv;
2116: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
2117: K = N; /* full rank */
2118: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
2119: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2121: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2122: Alpha = 1.0;
2123: ldb = lda;
2124: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2125: /* Ainv is Q, overwritten with inverse */
2127: if (debug) { /* Check that pseudo-inverse worked */
2128: PetscScalar Beta = 0.0;
2129: PetscBLASInt ldc;
2130: K = N;
2131: ldc = N;
2132: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2133: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2134: PetscFree(Aback);
2135: }
2136: return(0);
2137: }
2141: /* Overwrites A. Can handle degenerate problems and m<n. */
2142: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2143: {
2144: PetscBool debug = PETSC_FALSE;
2145: PetscScalar *Brhs, *Aback;
2146: PetscScalar *tmpwork;
2147: PetscReal rcond;
2148: #if defined (PETSC_USE_COMPLEX)
2149: PetscInt rworkSize;
2150: PetscReal *rwork;
2151: #endif
2152: PetscInt i, j, maxmn;
2153: PetscBLASInt M, N, lda, ldb, ldwork;
2154: PetscBLASInt nrhs, irank, info;
2158: if (debug) {
2159: PetscMalloc1(m*n,&Aback);
2160: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2161: }
2163: /* initialize to identity */
2164: tmpwork = Ainv;
2165: Brhs = work;
2166: maxmn = PetscMax(m,n);
2167: for (j=0; j<maxmn; j++) {
2168: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2169: }
2171: PetscBLASIntCast(m,&M);
2172: PetscBLASIntCast(n,&N);
2173: PetscBLASIntCast(mstride,&lda);
2174: PetscBLASIntCast(maxmn,&ldb);
2175: PetscBLASIntCast(worksize,&ldwork);
2176: rcond = -1;
2177: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2178: nrhs = M;
2179: #if defined(PETSC_USE_COMPLEX)
2180: rworkSize = 5 * PetscMin(M,N);
2181: PetscMalloc1(rworkSize,&rwork);
2182: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2183: PetscFPTrapPop();
2184: PetscFree(rwork);
2185: #else
2186: nrhs = M;
2187: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2188: PetscFPTrapPop();
2189: #endif
2190: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2191: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2192: if (irank < PetscMin(M,N)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2194: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2195: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2196: for (i=0; i<n; i++) {
2197: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2198: }
2199: return(0);
2200: }
2202: #if 0
2205: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2206: {
2207: PetscReal grad[2] = {0, 0};
2208: const PetscInt *faces;
2209: PetscInt numFaces, f;
2210: PetscErrorCode ierr;
2213: DMPlexGetConeSize(dm, cell, &numFaces);
2214: DMPlexGetCone(dm, cell, &faces);
2215: for (f = 0; f < numFaces; ++f) {
2216: const PetscInt *fcells;
2217: const CellGeom *cg1;
2218: const FaceGeom *fg;
2220: DMPlexGetSupport(dm, faces[f], &fcells);
2221: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2222: for (i = 0; i < 2; ++i) {
2223: PetscScalar du;
2225: if (fcells[i] == c) continue;
2226: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2227: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2228: grad[0] += fg->grad[!i][0] * du;
2229: grad[1] += fg->grad[!i][1] * du;
2230: }
2231: }
2232: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2233: return(0);
2234: }
2235: #endif
2239: /*
2240: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2242: Input Parameters:
2243: + fvm - The PetscFV object
2244: . numFaces - The number of cell faces which are not constrained
2245: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2247: Level: developer
2249: .seealso: PetscFVCreate()
2250: */
2251: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2252: {
2253: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2254: const PetscBool useSVD = PETSC_TRUE;
2255: const PetscInt maxFaces = ls->maxFaces;
2256: PetscInt dim, f, d;
2257: PetscErrorCode ierr;
2260: if (numFaces > maxFaces) {
2261: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2262: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2263: }
2264: PetscFVGetSpatialDimension(fvm, &dim);
2265: for (f = 0; f < numFaces; ++f) {
2266: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2267: }
2268: /* Overwrites B with garbage, returns Binv in row-major format */
2269: if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2270: else {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2271: for (f = 0; f < numFaces; ++f) {
2272: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2273: }
2274: return(0);
2275: }
2279: /*
2280: neighborVol[f*2+0] contains the left geom
2281: neighborVol[f*2+1] contains the right geom
2282: */
2283: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2284: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2285: {
2286: void (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
2287: void *rctx;
2288: PetscScalar *flux = fvm->fluxWork;
2289: PetscInt dim, pdim, Nc, totDim, off, f, d;
2293: PetscDSGetTotalComponents(prob, &Nc);
2294: PetscDSGetTotalDimension(prob, &totDim);
2295: PetscDSGetFieldOffset(prob, field, &off);
2296: PetscDSGetRiemannSolver(prob, field, &riemann);
2297: PetscDSGetContext(prob, field, &rctx);
2298: PetscFVGetSpatialDimension(fvm, &dim);
2299: PetscFVGetNumComponents(fvm, &pdim);
2300: for (f = 0; f < Nf; ++f) {
2301: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
2302: for (d = 0; d < pdim; ++d) {
2303: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2304: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2305: }
2306: }
2307: return(0);
2308: }
2312: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2313: {
2314: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2315: PetscInt dim, m, n, nrhs, minwork;
2316: PetscErrorCode ierr;
2320: PetscFVGetSpatialDimension(fvm, &dim);
2321: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2322: ls->maxFaces = maxFaces;
2323: m = ls->maxFaces;
2324: n = dim;
2325: nrhs = ls->maxFaces;
2326: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2327: ls->workSize = 5*minwork; /* We can afford to be extra generous */
2328: PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2329: return(0);
2330: }
2334: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2335: {
2337: fvm->ops->setfromoptions = NULL;
2338: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2339: fvm->ops->view = PetscFVView_LeastSquares;
2340: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2341: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2342: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2343: return(0);
2344: }
2346: /*MC
2347: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2349: Level: intermediate
2351: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2352: M*/
2356: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2357: {
2358: PetscFV_LeastSquares *ls;
2359: PetscErrorCode ierr;
2363: PetscNewLog(fvm, &ls);
2364: fvm->data = ls;
2366: ls->maxFaces = -1;
2367: ls->workSize = -1;
2368: ls->B = NULL;
2369: ls->Binv = NULL;
2370: ls->tau = NULL;
2371: ls->work = NULL;
2373: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2374: PetscFVInitialize_LeastSquares(fvm);
2375: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2376: return(0);
2377: }
2381: /*@
2382: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2384: Not collective
2386: Input parameters:
2387: + fvm - The PetscFV object
2388: - maxFaces - The maximum number of cell faces
2390: Level: intermediate
2392: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2393: @*/
2394: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2395: {
2400: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2401: return(0);
2402: }