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
24: Input Parameters:
25: + name - The name of a new user-defined creation routine
26: - create_func - The creation routine itself
28: Notes:
29: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
31: Sample usage:
32: .vb
33: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
34: .ve
36: Then, your PetscLimiter type can be chosen with the procedural interface via
37: .vb
38: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
39: PetscLimiterSetType(PetscLimiter, "my_lim");
40: .ve
41: or at runtime via the option
42: .vb
43: -petsclimiter_type my_lim
44: .ve
46: Level: advanced
48: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
50: @*/
51: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
52: {
53: PetscFunctionListAdd(&PetscLimiterList, sname, function);
54: return 0;
55: }
57: /*@C
58: PetscLimiterSetType - Builds a particular PetscLimiter
60: Collective on lim
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: PetscLimiterGetType(), PetscLimiterCreate()
72: @*/
73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
74: {
75: PetscErrorCode (*r)(PetscLimiter);
76: PetscBool match;
79: PetscObjectTypeCompare((PetscObject) lim, name, &match);
80: if (match) return 0;
82: PetscLimiterRegisterAll();
83: PetscFunctionListFind(PetscLimiterList, name, &r);
86: if (lim->ops->destroy) {
87: (*lim->ops->destroy)(lim);
88: lim->ops->destroy = NULL;
89: }
90: (*r)(lim);
91: PetscObjectChangeTypeName((PetscObject) lim, name);
92: return 0;
93: }
95: /*@C
96: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
98: Not Collective
100: Input Parameter:
101: . lim - The PetscLimiter
103: Output Parameter:
104: . name - The PetscLimiter type name
106: Level: intermediate
108: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
114: PetscLimiterRegisterAll();
115: *name = ((PetscObject) lim)->type_name;
116: return 0;
117: }
119: /*@C
120: PetscLimiterViewFromOptions - View from Options
122: Collective on PetscLimiter
124: Input Parameters:
125: + A - the PetscLimiter object to view
126: . obj - Optional object
127: - name - command line option
129: Level: intermediate
130: .seealso: PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
131: @*/
132: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
133: {
135: PetscObjectViewFromOptions((PetscObject)A,obj,name);
136: return 0;
137: }
139: /*@C
140: PetscLimiterView - Views a PetscLimiter
142: Collective on lim
144: Input Parameters:
145: + lim - the PetscLimiter object to view
146: - v - the viewer
148: Level: beginner
150: .seealso: PetscLimiterDestroy()
151: @*/
152: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
153: {
155: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);
156: if (lim->ops->view) (*lim->ops->view)(lim, v);
157: return 0;
158: }
160: /*@
161: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
163: Collective on lim
165: Input Parameter:
166: . lim - the PetscLimiter object to set options for
168: Level: intermediate
170: .seealso: PetscLimiterView()
171: @*/
172: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
173: {
174: const char *defaultType;
175: char name[256];
176: PetscBool flg;
180: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
181: else defaultType = ((PetscObject) lim)->type_name;
182: PetscLimiterRegisterAll();
184: PetscObjectOptionsBegin((PetscObject) lim);
185: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
186: if (flg) {
187: PetscLimiterSetType(lim, name);
188: } else if (!((PetscObject) lim)->type_name) {
189: PetscLimiterSetType(lim, defaultType);
190: }
191: if (lim->ops->setfromoptions) (*lim->ops->setfromoptions)(lim);
192: /* process any options handlers added with PetscObjectAddOptionsHandler() */
193: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
194: PetscOptionsEnd();
195: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
196: return 0;
197: }
199: /*@C
200: PetscLimiterSetUp - Construct data structures for the PetscLimiter
202: Collective on lim
204: Input Parameter:
205: . lim - the PetscLimiter object to setup
207: Level: intermediate
209: .seealso: PetscLimiterView(), PetscLimiterDestroy()
210: @*/
211: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
212: {
214: if (lim->ops->setup) (*lim->ops->setup)(lim);
215: return 0;
216: }
218: /*@
219: PetscLimiterDestroy - Destroys a PetscLimiter object
221: Collective on lim
223: Input Parameter:
224: . lim - the PetscLimiter object to destroy
226: Level: beginner
228: .seealso: PetscLimiterView()
229: @*/
230: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
231: {
232: if (!*lim) return 0;
235: if (--((PetscObject)(*lim))->refct > 0) {*lim = NULL; return 0;}
236: ((PetscObject) (*lim))->refct = 0;
238: if ((*lim)->ops->destroy) (*(*lim)->ops->destroy)(*lim);
239: PetscHeaderDestroy(lim);
240: return 0;
241: }
243: /*@
244: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
246: Collective
248: Input Parameter:
249: . comm - The communicator for the PetscLimiter object
251: Output Parameter:
252: . lim - The PetscLimiter object
254: Level: beginner
256: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
257: @*/
258: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
259: {
260: PetscLimiter l;
263: PetscCitationsRegister(LimiterCitation,&Limitercite);
264: *lim = NULL;
265: PetscFVInitializePackage();
267: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
269: *lim = l;
270: return 0;
271: }
273: /*@
274: PetscLimiterLimit - Limit the flux
276: Input Parameters:
277: + lim - The PetscLimiter
278: - flim - The input field
280: Output Parameter:
281: . phi - The limited field
283: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
284: $ The classical flux-limited formulation is psi(r) where
285: $
286: $ r = (u[0] - u[-1]) / (u[1] - u[0])
287: $
288: $ The second order TVD region is bounded by
289: $
290: $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
291: $
292: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
293: $ phi(r)(r+1)/2 in which the reconstructed interface values are
294: $
295: $ u(v) = u[0] + phi(r) (grad u)[0] v
296: $
297: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
298: $
299: $ 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))
300: $
301: $ For a nicer symmetric formulation, rewrite in terms of
302: $
303: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
304: $
305: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
306: $
307: $ phi(r) = phi(1/r)
308: $
309: $ becomes
310: $
311: $ w(f) = w(1-f).
312: $
313: $ The limiters below implement this final form w(f). The reference methods are
314: $
315: $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
317: Level: beginner
319: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
320: @*/
321: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
322: {
325: (*lim->ops->limit)(lim, flim, phi);
326: return 0;
327: }
329: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
330: {
331: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
333: PetscFree(l);
334: return 0;
335: }
337: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
338: {
339: PetscViewerFormat format;
341: PetscViewerGetFormat(viewer, &format);
342: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
343: return 0;
344: }
346: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
347: {
348: PetscBool iascii;
352: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
353: if (iascii) PetscLimiterView_Sin_Ascii(lim, viewer);
354: return 0;
355: }
357: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
358: {
359: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
360: return 0;
361: }
363: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
364: {
365: lim->ops->view = PetscLimiterView_Sin;
366: lim->ops->destroy = PetscLimiterDestroy_Sin;
367: lim->ops->limit = PetscLimiterLimit_Sin;
368: return 0;
369: }
371: /*MC
372: PETSCLIMITERSIN = "sin" - A PetscLimiter object
374: Level: intermediate
376: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
377: M*/
379: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
380: {
381: PetscLimiter_Sin *l;
384: PetscNewLog(lim, &l);
385: lim->data = l;
387: PetscLimiterInitialize_Sin(lim);
388: return 0;
389: }
391: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
392: {
393: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
395: PetscFree(l);
396: return 0;
397: }
399: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
400: {
401: PetscViewerFormat format;
403: PetscViewerGetFormat(viewer, &format);
404: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
405: return 0;
406: }
408: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
409: {
410: PetscBool iascii;
414: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
415: if (iascii) PetscLimiterView_Zero_Ascii(lim, viewer);
416: return 0;
417: }
419: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
420: {
421: *phi = 0.0;
422: return 0;
423: }
425: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
426: {
427: lim->ops->view = PetscLimiterView_Zero;
428: lim->ops->destroy = PetscLimiterDestroy_Zero;
429: lim->ops->limit = PetscLimiterLimit_Zero;
430: return 0;
431: }
433: /*MC
434: PETSCLIMITERZERO = "zero" - A PetscLimiter object
436: Level: intermediate
438: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
439: M*/
441: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
442: {
443: PetscLimiter_Zero *l;
446: PetscNewLog(lim, &l);
447: lim->data = l;
449: PetscLimiterInitialize_Zero(lim);
450: return 0;
451: }
453: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
454: {
455: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
457: PetscFree(l);
458: return 0;
459: }
461: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
462: {
463: PetscViewerFormat format;
465: PetscViewerGetFormat(viewer, &format);
466: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
467: return 0;
468: }
470: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
471: {
472: PetscBool iascii;
476: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
477: if (iascii) PetscLimiterView_None_Ascii(lim, viewer);
478: return 0;
479: }
481: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
482: {
483: *phi = 1.0;
484: return 0;
485: }
487: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
488: {
489: lim->ops->view = PetscLimiterView_None;
490: lim->ops->destroy = PetscLimiterDestroy_None;
491: lim->ops->limit = PetscLimiterLimit_None;
492: return 0;
493: }
495: /*MC
496: PETSCLIMITERNONE = "none" - A PetscLimiter object
498: Level: intermediate
500: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
501: M*/
503: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
504: {
505: PetscLimiter_None *l;
508: PetscNewLog(lim, &l);
509: lim->data = l;
511: PetscLimiterInitialize_None(lim);
512: return 0;
513: }
515: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
516: {
517: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
519: PetscFree(l);
520: return 0;
521: }
523: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
524: {
525: PetscViewerFormat format;
527: PetscViewerGetFormat(viewer, &format);
528: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
529: return 0;
530: }
532: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
533: {
534: PetscBool iascii;
538: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
539: if (iascii) PetscLimiterView_Minmod_Ascii(lim, viewer);
540: return 0;
541: }
543: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
544: {
545: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
546: return 0;
547: }
549: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
550: {
551: lim->ops->view = PetscLimiterView_Minmod;
552: lim->ops->destroy = PetscLimiterDestroy_Minmod;
553: lim->ops->limit = PetscLimiterLimit_Minmod;
554: return 0;
555: }
557: /*MC
558: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
560: Level: intermediate
562: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
563: M*/
565: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
566: {
567: PetscLimiter_Minmod *l;
570: PetscNewLog(lim, &l);
571: lim->data = l;
573: PetscLimiterInitialize_Minmod(lim);
574: return 0;
575: }
577: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
578: {
579: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
581: PetscFree(l);
582: return 0;
583: }
585: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
586: {
587: PetscViewerFormat format;
589: PetscViewerGetFormat(viewer, &format);
590: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
591: return 0;
592: }
594: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
595: {
596: PetscBool iascii;
600: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
601: if (iascii) PetscLimiterView_VanLeer_Ascii(lim, viewer);
602: return 0;
603: }
605: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
606: {
607: *phi = PetscMax(0, 4*f*(1-f));
608: return 0;
609: }
611: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
612: {
613: lim->ops->view = PetscLimiterView_VanLeer;
614: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
615: lim->ops->limit = PetscLimiterLimit_VanLeer;
616: return 0;
617: }
619: /*MC
620: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
622: Level: intermediate
624: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
625: M*/
627: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
628: {
629: PetscLimiter_VanLeer *l;
632: PetscNewLog(lim, &l);
633: lim->data = l;
635: PetscLimiterInitialize_VanLeer(lim);
636: return 0;
637: }
639: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
640: {
641: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
643: PetscFree(l);
644: return 0;
645: }
647: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
648: {
649: PetscViewerFormat format;
651: PetscViewerGetFormat(viewer, &format);
652: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
653: return 0;
654: }
656: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
657: {
658: PetscBool iascii;
662: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
663: if (iascii) PetscLimiterView_VanAlbada_Ascii(lim, viewer);
664: return 0;
665: }
667: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
668: {
669: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
670: return 0;
671: }
673: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
674: {
675: lim->ops->view = PetscLimiterView_VanAlbada;
676: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
677: lim->ops->limit = PetscLimiterLimit_VanAlbada;
678: return 0;
679: }
681: /*MC
682: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
684: Level: intermediate
686: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
687: M*/
689: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
690: {
691: PetscLimiter_VanAlbada *l;
694: PetscNewLog(lim, &l);
695: lim->data = l;
697: PetscLimiterInitialize_VanAlbada(lim);
698: return 0;
699: }
701: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
702: {
703: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
705: PetscFree(l);
706: return 0;
707: }
709: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
710: {
711: PetscViewerFormat format;
713: PetscViewerGetFormat(viewer, &format);
714: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
715: return 0;
716: }
718: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
719: {
720: PetscBool iascii;
724: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
725: if (iascii) PetscLimiterView_Superbee_Ascii(lim, viewer);
726: return 0;
727: }
729: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
730: {
731: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
732: return 0;
733: }
735: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
736: {
737: lim->ops->view = PetscLimiterView_Superbee;
738: lim->ops->destroy = PetscLimiterDestroy_Superbee;
739: lim->ops->limit = PetscLimiterLimit_Superbee;
740: return 0;
741: }
743: /*MC
744: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
746: Level: intermediate
748: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
749: M*/
751: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
752: {
753: PetscLimiter_Superbee *l;
756: PetscNewLog(lim, &l);
757: lim->data = l;
759: PetscLimiterInitialize_Superbee(lim);
760: return 0;
761: }
763: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
764: {
765: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
767: PetscFree(l);
768: return 0;
769: }
771: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
772: {
773: PetscViewerFormat format;
775: PetscViewerGetFormat(viewer, &format);
776: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
777: return 0;
778: }
780: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
781: {
782: PetscBool iascii;
786: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
787: if (iascii) PetscLimiterView_MC_Ascii(lim, viewer);
788: return 0;
789: }
791: /* aka Barth-Jespersen */
792: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
793: {
794: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
795: return 0;
796: }
798: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
799: {
800: lim->ops->view = PetscLimiterView_MC;
801: lim->ops->destroy = PetscLimiterDestroy_MC;
802: lim->ops->limit = PetscLimiterLimit_MC;
803: return 0;
804: }
806: /*MC
807: PETSCLIMITERMC = "mc" - A PetscLimiter object
809: Level: intermediate
811: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
812: M*/
814: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
815: {
816: PetscLimiter_MC *l;
819: PetscNewLog(lim, &l);
820: lim->data = l;
822: PetscLimiterInitialize_MC(lim);
823: return 0;
824: }
826: PetscClassId PETSCFV_CLASSID = 0;
828: PetscFunctionList PetscFVList = NULL;
829: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
831: /*@C
832: PetscFVRegister - Adds a new PetscFV implementation
834: Not Collective
836: Input Parameters:
837: + name - The name of a new user-defined creation routine
838: - create_func - The creation routine itself
840: Notes:
841: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
843: Sample usage:
844: .vb
845: PetscFVRegister("my_fv", MyPetscFVCreate);
846: .ve
848: Then, your PetscFV type can be chosen with the procedural interface via
849: .vb
850: PetscFVCreate(MPI_Comm, PetscFV *);
851: PetscFVSetType(PetscFV, "my_fv");
852: .ve
853: or at runtime via the option
854: .vb
855: -petscfv_type my_fv
856: .ve
858: Level: advanced
860: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
862: @*/
863: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
864: {
865: PetscFunctionListAdd(&PetscFVList, sname, function);
866: return 0;
867: }
869: /*@C
870: PetscFVSetType - Builds a particular PetscFV
872: Collective on fvm
874: Input Parameters:
875: + fvm - The PetscFV object
876: - name - The kind of FVM space
878: Options Database Key:
879: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
881: Level: intermediate
883: .seealso: PetscFVGetType(), PetscFVCreate()
884: @*/
885: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
886: {
887: PetscErrorCode (*r)(PetscFV);
888: PetscBool match;
891: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
892: if (match) return 0;
894: PetscFVRegisterAll();
895: PetscFunctionListFind(PetscFVList, name, &r);
898: if (fvm->ops->destroy) {
899: (*fvm->ops->destroy)(fvm);
900: fvm->ops->destroy = NULL;
901: }
902: (*r)(fvm);
903: PetscObjectChangeTypeName((PetscObject) fvm, name);
904: return 0;
905: }
907: /*@C
908: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
910: Not Collective
912: Input Parameter:
913: . fvm - The PetscFV
915: Output Parameter:
916: . name - The PetscFV type name
918: Level: intermediate
920: .seealso: PetscFVSetType(), PetscFVCreate()
921: @*/
922: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
923: {
926: PetscFVRegisterAll();
927: *name = ((PetscObject) fvm)->type_name;
928: return 0;
929: }
931: /*@C
932: PetscFVViewFromOptions - View from Options
934: Collective on PetscFV
936: Input Parameters:
937: + A - the PetscFV object
938: . obj - Optional object
939: - name - command line option
941: Level: intermediate
942: .seealso: PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
943: @*/
944: PetscErrorCode PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
945: {
947: PetscObjectViewFromOptions((PetscObject)A,obj,name);
948: return 0;
949: }
951: /*@C
952: PetscFVView - Views a PetscFV
954: Collective on fvm
956: Input Parameters:
957: + fvm - the PetscFV object to view
958: - v - the viewer
960: Level: beginner
962: .seealso: PetscFVDestroy()
963: @*/
964: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
965: {
967: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);
968: if (fvm->ops->view) (*fvm->ops->view)(fvm, v);
969: return 0;
970: }
972: /*@
973: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
975: Collective on fvm
977: Input Parameter:
978: . fvm - the PetscFV object to set options for
980: Options Database Key:
981: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
983: Level: intermediate
985: .seealso: PetscFVView()
986: @*/
987: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
988: {
989: const char *defaultType;
990: char name[256];
991: PetscBool flg;
995: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
996: else defaultType = ((PetscObject) fvm)->type_name;
997: PetscFVRegisterAll();
999: PetscObjectOptionsBegin((PetscObject) fvm);
1000: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1001: if (flg) {
1002: PetscFVSetType(fvm, name);
1003: } else if (!((PetscObject) fvm)->type_name) {
1004: PetscFVSetType(fvm, defaultType);
1006: }
1007: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1008: if (fvm->ops->setfromoptions) (*fvm->ops->setfromoptions)(fvm);
1009: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1010: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1011: PetscLimiterSetFromOptions(fvm->limiter);
1012: PetscOptionsEnd();
1013: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1014: return 0;
1015: }
1017: /*@
1018: PetscFVSetUp - Construct data structures for the PetscFV
1020: Collective on fvm
1022: Input Parameter:
1023: . fvm - the PetscFV object to setup
1025: Level: intermediate
1027: .seealso: PetscFVView(), PetscFVDestroy()
1028: @*/
1029: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1030: {
1032: PetscLimiterSetUp(fvm->limiter);
1033: if (fvm->ops->setup) (*fvm->ops->setup)(fvm);
1034: return 0;
1035: }
1037: /*@
1038: PetscFVDestroy - Destroys a PetscFV object
1040: Collective on fvm
1042: Input Parameter:
1043: . fvm - the PetscFV object to destroy
1045: Level: beginner
1047: .seealso: PetscFVView()
1048: @*/
1049: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1050: {
1051: PetscInt i;
1053: if (!*fvm) return 0;
1056: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; return 0;}
1057: ((PetscObject) (*fvm))->refct = 0;
1059: for (i = 0; i < (*fvm)->numComponents; i++) {
1060: PetscFree((*fvm)->componentNames[i]);
1061: }
1062: PetscFree((*fvm)->componentNames);
1063: PetscLimiterDestroy(&(*fvm)->limiter);
1064: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1065: PetscFree((*fvm)->fluxWork);
1066: PetscQuadratureDestroy(&(*fvm)->quadrature);
1067: PetscTabulationDestroy(&(*fvm)->T);
1069: if ((*fvm)->ops->destroy) (*(*fvm)->ops->destroy)(*fvm);
1070: PetscHeaderDestroy(fvm);
1071: return 0;
1072: }
1074: /*@
1075: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1077: Collective
1079: Input Parameter:
1080: . comm - The communicator for the PetscFV object
1082: Output Parameter:
1083: . fvm - The PetscFV object
1085: Level: beginner
1087: .seealso: PetscFVSetType(), PETSCFVUPWIND
1088: @*/
1089: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1090: {
1091: PetscFV f;
1094: *fvm = NULL;
1095: PetscFVInitializePackage();
1097: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1098: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1100: PetscLimiterCreate(comm, &f->limiter);
1101: f->numComponents = 1;
1102: f->dim = 0;
1103: f->computeGradients = PETSC_FALSE;
1104: f->fluxWork = NULL;
1105: PetscCalloc1(f->numComponents,&f->componentNames);
1107: *fvm = f;
1108: return 0;
1109: }
1111: /*@
1112: PetscFVSetLimiter - Set the limiter object
1114: Logically collective on fvm
1116: Input Parameters:
1117: + fvm - the PetscFV object
1118: - lim - The PetscLimiter
1120: Level: intermediate
1122: .seealso: PetscFVGetLimiter()
1123: @*/
1124: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1125: {
1128: PetscLimiterDestroy(&fvm->limiter);
1129: PetscObjectReference((PetscObject) lim);
1130: fvm->limiter = lim;
1131: return 0;
1132: }
1134: /*@
1135: PetscFVGetLimiter - Get the limiter object
1137: Not collective
1139: Input Parameter:
1140: . fvm - the PetscFV object
1142: Output Parameter:
1143: . lim - The PetscLimiter
1145: Level: intermediate
1147: .seealso: PetscFVSetLimiter()
1148: @*/
1149: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1150: {
1153: *lim = fvm->limiter;
1154: return 0;
1155: }
1157: /*@
1158: PetscFVSetNumComponents - Set the number of field components
1160: Logically collective on fvm
1162: Input Parameters:
1163: + fvm - the PetscFV object
1164: - comp - The number of components
1166: Level: intermediate
1168: .seealso: PetscFVGetNumComponents()
1169: @*/
1170: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1171: {
1173: if (fvm->numComponents != comp) {
1174: PetscInt i;
1176: for (i = 0; i < fvm->numComponents; i++) {
1177: PetscFree(fvm->componentNames[i]);
1178: }
1179: PetscFree(fvm->componentNames);
1180: PetscCalloc1(comp,&fvm->componentNames);
1181: }
1182: fvm->numComponents = comp;
1183: PetscFree(fvm->fluxWork);
1184: PetscMalloc1(comp, &fvm->fluxWork);
1185: return 0;
1186: }
1188: /*@
1189: PetscFVGetNumComponents - Get the number of field components
1191: Not collective
1193: Input Parameter:
1194: . fvm - the PetscFV object
1196: Output Parameter:
1197: , comp - The number of components
1199: Level: intermediate
1201: .seealso: PetscFVSetNumComponents()
1202: @*/
1203: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1204: {
1207: *comp = fvm->numComponents;
1208: return 0;
1209: }
1211: /*@C
1212: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1214: Logically collective on fvm
1215: Input Parameters:
1216: + fvm - the PetscFV object
1217: . comp - the component number
1218: - name - the component name
1220: Level: intermediate
1222: .seealso: PetscFVGetComponentName()
1223: @*/
1224: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1225: {
1226: PetscFree(fvm->componentNames[comp]);
1227: PetscStrallocpy(name,&fvm->componentNames[comp]);
1228: return 0;
1229: }
1231: /*@C
1232: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1234: Logically collective on fvm
1235: Input Parameters:
1236: + fvm - the PetscFV object
1237: - comp - the component number
1239: Output Parameter:
1240: . name - the component name
1242: Level: intermediate
1244: .seealso: PetscFVSetComponentName()
1245: @*/
1246: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1247: {
1248: *name = fvm->componentNames[comp];
1249: return 0;
1250: }
1252: /*@
1253: PetscFVSetSpatialDimension - Set the spatial dimension
1255: Logically collective on fvm
1257: Input Parameters:
1258: + fvm - the PetscFV object
1259: - dim - The spatial dimension
1261: Level: intermediate
1263: .seealso: PetscFVGetSpatialDimension()
1264: @*/
1265: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1266: {
1268: fvm->dim = dim;
1269: return 0;
1270: }
1272: /*@
1273: PetscFVGetSpatialDimension - Get the spatial dimension
1275: Logically collective on fvm
1277: Input Parameter:
1278: . fvm - the PetscFV object
1280: Output Parameter:
1281: . dim - The spatial dimension
1283: Level: intermediate
1285: .seealso: PetscFVSetSpatialDimension()
1286: @*/
1287: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1288: {
1291: *dim = fvm->dim;
1292: return 0;
1293: }
1295: /*@
1296: PetscFVSetComputeGradients - Toggle computation of cell gradients
1298: Logically collective on fvm
1300: Input Parameters:
1301: + fvm - the PetscFV object
1302: - computeGradients - Flag to compute cell gradients
1304: Level: intermediate
1306: .seealso: PetscFVGetComputeGradients()
1307: @*/
1308: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1309: {
1311: fvm->computeGradients = computeGradients;
1312: return 0;
1313: }
1315: /*@
1316: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1318: Not collective
1320: Input Parameter:
1321: . fvm - the PetscFV object
1323: Output Parameter:
1324: . computeGradients - Flag to compute cell gradients
1326: Level: intermediate
1328: .seealso: PetscFVSetComputeGradients()
1329: @*/
1330: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1331: {
1334: *computeGradients = fvm->computeGradients;
1335: return 0;
1336: }
1338: /*@
1339: PetscFVSetQuadrature - Set the quadrature object
1341: Logically collective on fvm
1343: Input Parameters:
1344: + fvm - the PetscFV object
1345: - q - The PetscQuadrature
1347: Level: intermediate
1349: .seealso: PetscFVGetQuadrature()
1350: @*/
1351: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1352: {
1354: PetscQuadratureDestroy(&fvm->quadrature);
1355: PetscObjectReference((PetscObject) q);
1356: fvm->quadrature = q;
1357: return 0;
1358: }
1360: /*@
1361: PetscFVGetQuadrature - Get the quadrature object
1363: Not collective
1365: Input Parameter:
1366: . fvm - the PetscFV object
1368: Output Parameter:
1369: . lim - The PetscQuadrature
1371: Level: intermediate
1373: .seealso: PetscFVSetQuadrature()
1374: @*/
1375: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1376: {
1379: if (!fvm->quadrature) {
1380: /* Create default 1-point quadrature */
1381: PetscReal *points, *weights;
1383: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1384: PetscCalloc1(fvm->dim, &points);
1385: PetscMalloc1(1, &weights);
1386: weights[0] = 1.0;
1387: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1388: }
1389: *q = fvm->quadrature;
1390: return 0;
1391: }
1393: /*@
1394: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1396: Not collective
1398: Input Parameter:
1399: . fvm - The PetscFV object
1401: Output Parameter:
1402: . sp - The PetscDualSpace object
1404: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1406: Level: intermediate
1408: .seealso: PetscFVCreate()
1409: @*/
1410: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1411: {
1414: if (!fvm->dualSpace) {
1415: DM K;
1416: PetscInt dim, Nc, c;
1418: PetscFVGetSpatialDimension(fvm, &dim);
1419: PetscFVGetNumComponents(fvm, &Nc);
1420: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1421: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1422: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K);
1423: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1424: PetscDualSpaceSetDM(fvm->dualSpace, K);
1425: DMDestroy(&K);
1426: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1427: /* Should we be using PetscFVGetQuadrature() here? */
1428: for (c = 0; c < Nc; ++c) {
1429: PetscQuadrature qc;
1430: PetscReal *points, *weights;
1432: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1433: PetscCalloc1(dim, &points);
1434: PetscCalloc1(Nc, &weights);
1435: weights[c] = 1.0;
1436: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1437: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1438: PetscQuadratureDestroy(&qc);
1439: }
1440: PetscDualSpaceSetUp(fvm->dualSpace);
1441: }
1442: *sp = fvm->dualSpace;
1443: return 0;
1444: }
1446: /*@
1447: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1449: Not collective
1451: Input Parameters:
1452: + fvm - The PetscFV object
1453: - sp - The PetscDualSpace object
1455: Level: intermediate
1457: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1459: .seealso: PetscFVCreate()
1460: @*/
1461: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1462: {
1465: PetscDualSpaceDestroy(&fvm->dualSpace);
1466: fvm->dualSpace = sp;
1467: PetscObjectReference((PetscObject) fvm->dualSpace);
1468: return 0;
1469: }
1471: /*@C
1472: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1474: Not collective
1476: Input Parameter:
1477: . fvm - The PetscFV object
1479: Output Parameter:
1480: . T - The basis function values and derivatives at quadrature points
1482: Note:
1483: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1484: $ 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
1485: $ 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
1487: Level: intermediate
1489: .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1490: @*/
1491: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1492: {
1493: PetscInt npoints;
1494: const PetscReal *points;
1498: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1499: if (!fvm->T) PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);
1500: *T = fvm->T;
1501: return 0;
1502: }
1504: /*@C
1505: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1507: Not collective
1509: Input Parameters:
1510: + fvm - The PetscFV object
1511: . nrepl - The number of replicas
1512: . npoints - The number of tabulation points in a replica
1513: . points - The tabulation point coordinates
1514: - K - The order of derivative to tabulate
1516: Output Parameter:
1517: . T - The basis function values and derivative at tabulation points
1519: Note:
1520: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1521: $ 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
1522: $ 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
1524: Level: intermediate
1526: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1527: @*/
1528: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1529: {
1530: PetscInt pdim = 1; /* Dimension of approximation space P */
1531: PetscInt cdim; /* Spatial dimension */
1532: PetscInt Nc; /* Field components */
1533: PetscInt k, p, d, c, e;
1535: if (!npoints || K < 0) {
1536: *T = NULL;
1537: return 0;
1538: }
1542: PetscFVGetSpatialDimension(fvm, &cdim);
1543: PetscFVGetNumComponents(fvm, &Nc);
1544: PetscMalloc1(1, T);
1545: (*T)->K = !cdim ? 0 : K;
1546: (*T)->Nr = nrepl;
1547: (*T)->Np = npoints;
1548: (*T)->Nb = pdim;
1549: (*T)->Nc = Nc;
1550: (*T)->cdim = cdim;
1551: PetscMalloc1((*T)->K+1, &(*T)->T);
1552: for (k = 0; k <= (*T)->K; ++k) {
1553: PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
1554: }
1555: if (K >= 0) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) (*T)->T[0][(p*pdim + d)*Nc + c] = 1.0;}
1556: if (K >= 1) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim; ++e) (*T)->T[1][((p*pdim + d)*Nc + c)*cdim + e] = 0.0;}
1557: if (K >= 2) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim*cdim; ++e) (*T)->T[2][((p*pdim + d)*Nc + c)*cdim*cdim + e] = 0.0;}
1558: return 0;
1559: }
1561: /*@C
1562: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1564: Input Parameters:
1565: + fvm - The PetscFV object
1566: . numFaces - The number of cell faces which are not constrained
1567: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1569: Level: advanced
1571: .seealso: PetscFVCreate()
1572: @*/
1573: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1574: {
1576: if (fvm->ops->computegradient) (*fvm->ops->computegradient)(fvm, numFaces, dx, grad);
1577: return 0;
1578: }
1580: /*@C
1581: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1583: Not collective
1585: Input Parameters:
1586: + fvm - The PetscFV object for the field being integrated
1587: . prob - The PetscDS specifing the discretizations and continuum functions
1588: . field - The field being integrated
1589: . Nf - The number of faces in the chunk
1590: . fgeom - The face geometry for each face in the chunk
1591: . neighborVol - The volume for each pair of cells in the chunk
1592: . uL - The state from the cell on the left
1593: - uR - The state from the cell on the right
1595: Output Parameters:
1596: + fluxL - the left fluxes for each face
1597: - fluxR - the right fluxes for each face
1599: Level: developer
1601: .seealso: PetscFVCreate()
1602: @*/
1603: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1604: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1605: {
1607: if (fvm->ops->integraterhsfunction) (*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1608: return 0;
1609: }
1611: /*@
1612: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1613: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1614: sparsity). It is also used to create an interpolation between regularly refined meshes.
1616: Input Parameter:
1617: . fv - The initial PetscFV
1619: Output Parameter:
1620: . fvRef - The refined PetscFV
1622: Level: advanced
1624: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1625: @*/
1626: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1627: {
1628: PetscDualSpace Q, Qref;
1629: DM K, Kref;
1630: PetscQuadrature q, qref;
1631: DMPolytopeType ct;
1632: DMPlexTransform tr;
1633: PetscReal *v0;
1634: PetscReal *jac, *invjac;
1635: PetscInt numComp, numSubelements, s;
1637: PetscFVGetDualSpace(fv, &Q);
1638: PetscFVGetQuadrature(fv, &q);
1639: PetscDualSpaceGetDM(Q, &K);
1640: /* Create dual space */
1641: PetscDualSpaceDuplicate(Q, &Qref);
1642: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1643: PetscDualSpaceSetDM(Qref, Kref);
1644: DMDestroy(&Kref);
1645: PetscDualSpaceSetUp(Qref);
1646: /* Create volume */
1647: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1648: PetscFVSetDualSpace(*fvRef, Qref);
1649: PetscFVGetNumComponents(fv, &numComp);
1650: PetscFVSetNumComponents(*fvRef, numComp);
1651: PetscFVSetUp(*fvRef);
1652: /* Create quadrature */
1653: DMPlexGetCellType(K, 0, &ct);
1654: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1655: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1656: DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1657: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1658: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1659: for (s = 0; s < numSubelements; ++s) {
1660: PetscQuadrature qs;
1661: const PetscReal *points, *weights;
1662: PetscReal *p, *w;
1663: PetscInt dim, Nc, npoints, np;
1665: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1666: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1667: np = npoints/numSubelements;
1668: PetscMalloc1(np*dim,&p);
1669: PetscMalloc1(np*Nc,&w);
1670: PetscArraycpy(p, &points[s*np*dim], np*dim);
1671: PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1672: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1673: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1674: PetscQuadratureDestroy(&qs);
1675: }
1676: PetscFVSetQuadrature(*fvRef, qref);
1677: DMPlexTransformDestroy(&tr);
1678: PetscQuadratureDestroy(&qref);
1679: PetscDualSpaceDestroy(&Qref);
1680: return 0;
1681: }
1683: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1684: {
1685: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1687: PetscFree(b);
1688: return 0;
1689: }
1691: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1692: {
1693: PetscInt Nc, c;
1694: PetscViewerFormat format;
1696: PetscFVGetNumComponents(fv, &Nc);
1697: PetscViewerGetFormat(viewer, &format);
1698: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1699: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1700: for (c = 0; c < Nc; c++) {
1701: if (fv->componentNames[c]) {
1702: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1703: }
1704: }
1705: return 0;
1706: }
1708: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1709: {
1710: PetscBool iascii;
1714: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1715: if (iascii) PetscFVView_Upwind_Ascii(fv, viewer);
1716: return 0;
1717: }
1719: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1720: {
1721: return 0;
1722: }
1724: /*
1725: neighborVol[f*2+0] contains the left geom
1726: neighborVol[f*2+1] contains the right geom
1727: */
1728: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1729: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1730: {
1731: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1732: void *rctx;
1733: PetscScalar *flux = fvm->fluxWork;
1734: const PetscScalar *constants;
1735: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1737: PetscDSGetTotalComponents(prob, &Nc);
1738: PetscDSGetTotalDimension(prob, &totDim);
1739: PetscDSGetFieldOffset(prob, field, &off);
1740: PetscDSGetRiemannSolver(prob, field, &riemann);
1741: PetscDSGetContext(prob, field, &rctx);
1742: PetscDSGetConstants(prob, &numConstants, &constants);
1743: PetscFVGetSpatialDimension(fvm, &dim);
1744: PetscFVGetNumComponents(fvm, &pdim);
1745: for (f = 0; f < Nf; ++f) {
1746: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1747: for (d = 0; d < pdim; ++d) {
1748: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1749: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1750: }
1751: }
1752: return 0;
1753: }
1755: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1756: {
1757: fvm->ops->setfromoptions = NULL;
1758: fvm->ops->setup = PetscFVSetUp_Upwind;
1759: fvm->ops->view = PetscFVView_Upwind;
1760: fvm->ops->destroy = PetscFVDestroy_Upwind;
1761: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1762: return 0;
1763: }
1765: /*MC
1766: PETSCFVUPWIND = "upwind" - A PetscFV object
1768: Level: intermediate
1770: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1771: M*/
1773: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1774: {
1775: PetscFV_Upwind *b;
1778: PetscNewLog(fvm,&b);
1779: fvm->data = b;
1781: PetscFVInitialize_Upwind(fvm);
1782: return 0;
1783: }
1785: #include <petscblaslapack.h>
1787: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1788: {
1789: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1791: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1792: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1793: PetscFree(ls);
1794: return 0;
1795: }
1797: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1798: {
1799: PetscInt Nc, c;
1800: PetscViewerFormat format;
1802: PetscFVGetNumComponents(fv, &Nc);
1803: PetscViewerGetFormat(viewer, &format);
1804: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1805: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1806: for (c = 0; c < Nc; c++) {
1807: if (fv->componentNames[c]) {
1808: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1809: }
1810: }
1811: return 0;
1812: }
1814: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1815: {
1816: PetscBool iascii;
1820: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1821: if (iascii) PetscFVView_LeastSquares_Ascii(fv, viewer);
1822: return 0;
1823: }
1825: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1826: {
1827: return 0;
1828: }
1830: /* Overwrites A. Can only handle full-rank problems with m>=n */
1831: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1832: {
1833: PetscBool debug = PETSC_FALSE;
1834: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1835: PetscScalar *R,*Q,*Aback,Alpha;
1837: if (debug) {
1838: PetscMalloc1(m*n,&Aback);
1839: PetscArraycpy(Aback,A,m*n);
1840: }
1842: PetscBLASIntCast(m,&M);
1843: PetscBLASIntCast(n,&N);
1844: PetscBLASIntCast(mstride,&lda);
1845: PetscBLASIntCast(worksize,&ldwork);
1846: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1847: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1848: PetscFPTrapPop();
1850: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1852: /* Extract an explicit representation of Q */
1853: Q = Ainv;
1854: PetscArraycpy(Q,A,mstride*n);
1855: K = N; /* full rank */
1856: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1859: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1860: Alpha = 1.0;
1861: ldb = lda;
1862: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1863: /* Ainv is Q, overwritten with inverse */
1865: if (debug) { /* Check that pseudo-inverse worked */
1866: PetscScalar Beta = 0.0;
1867: PetscBLASInt ldc;
1868: K = N;
1869: ldc = N;
1870: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1871: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1872: PetscFree(Aback);
1873: }
1874: return 0;
1875: }
1877: /* Overwrites A. Can handle degenerate problems and m<n. */
1878: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1879: {
1880: PetscBool debug = PETSC_FALSE;
1881: PetscScalar *Brhs, *Aback;
1882: PetscScalar *tmpwork;
1883: PetscReal rcond;
1884: #if defined (PETSC_USE_COMPLEX)
1885: PetscInt rworkSize;
1886: PetscReal *rwork;
1887: #endif
1888: PetscInt i, j, maxmn;
1889: PetscBLASInt M, N, lda, ldb, ldwork;
1890: PetscBLASInt nrhs, irank, info;
1892: if (debug) {
1893: PetscMalloc1(m*n,&Aback);
1894: PetscArraycpy(Aback,A,m*n);
1895: }
1897: /* initialize to identity */
1898: tmpwork = work;
1899: Brhs = Ainv;
1900: maxmn = PetscMax(m,n);
1901: for (j=0; j<maxmn; j++) {
1902: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1903: }
1905: PetscBLASIntCast(m,&M);
1906: PetscBLASIntCast(n,&N);
1907: PetscBLASIntCast(mstride,&lda);
1908: PetscBLASIntCast(maxmn,&ldb);
1909: PetscBLASIntCast(worksize,&ldwork);
1910: rcond = -1;
1911: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1912: nrhs = M;
1913: #if defined(PETSC_USE_COMPLEX)
1914: rworkSize = 5 * PetscMin(M,N);
1915: PetscMalloc1(rworkSize,&rwork);
1916: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info));
1917: PetscFPTrapPop();
1918: PetscFree(rwork);
1919: #else
1920: nrhs = M;
1921: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info));
1922: PetscFPTrapPop();
1923: #endif
1925: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1927: return 0;
1928: }
1930: #if 0
1931: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1932: {
1933: PetscReal grad[2] = {0, 0};
1934: const PetscInt *faces;
1935: PetscInt numFaces, f;
1937: DMPlexGetConeSize(dm, cell, &numFaces);
1938: DMPlexGetCone(dm, cell, &faces);
1939: for (f = 0; f < numFaces; ++f) {
1940: const PetscInt *fcells;
1941: const CellGeom *cg1;
1942: const FaceGeom *fg;
1944: DMPlexGetSupport(dm, faces[f], &fcells);
1945: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1946: for (i = 0; i < 2; ++i) {
1947: PetscScalar du;
1949: if (fcells[i] == c) continue;
1950: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1951: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1952: grad[0] += fg->grad[!i][0] * du;
1953: grad[1] += fg->grad[!i][1] * du;
1954: }
1955: }
1956: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1957: return 0;
1958: }
1959: #endif
1961: /*
1962: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1964: Input Parameters:
1965: + fvm - The PetscFV object
1966: . numFaces - The number of cell faces which are not constrained
1967: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1969: Level: developer
1971: .seealso: PetscFVCreate()
1972: */
1973: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1974: {
1975: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1976: const PetscBool useSVD = PETSC_TRUE;
1977: const PetscInt maxFaces = ls->maxFaces;
1978: PetscInt dim, f, d;
1980: if (numFaces > maxFaces) {
1982: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
1983: }
1984: PetscFVGetSpatialDimension(fvm, &dim);
1985: for (f = 0; f < numFaces; ++f) {
1986: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
1987: }
1988: /* Overwrites B with garbage, returns Binv in row-major format */
1989: if (useSVD) {
1990: PetscInt maxmn = PetscMax(numFaces, dim);
1991: PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1992: /* Binv shaped in column-major, coldim=maxmn.*/
1993: for (f = 0; f < numFaces; ++f) {
1994: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d + maxmn*f];
1995: }
1996: } else {
1997: PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1998: /* Binv shaped in row-major, rowdim=maxFaces.*/
1999: for (f = 0; f < numFaces; ++f) {
2000: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces + f];
2001: }
2002: }
2003: return 0;
2004: }
2006: /*
2007: neighborVol[f*2+0] contains the left geom
2008: neighborVol[f*2+1] contains the right geom
2009: */
2010: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2011: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2012: {
2013: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2014: void *rctx;
2015: PetscScalar *flux = fvm->fluxWork;
2016: const PetscScalar *constants;
2017: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2019: PetscDSGetTotalComponents(prob, &Nc);
2020: PetscDSGetTotalDimension(prob, &totDim);
2021: PetscDSGetFieldOffset(prob, field, &off);
2022: PetscDSGetRiemannSolver(prob, field, &riemann);
2023: PetscDSGetContext(prob, field, &rctx);
2024: PetscDSGetConstants(prob, &numConstants, &constants);
2025: PetscFVGetSpatialDimension(fvm, &dim);
2026: PetscFVGetNumComponents(fvm, &pdim);
2027: for (f = 0; f < Nf; ++f) {
2028: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2029: for (d = 0; d < pdim; ++d) {
2030: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2031: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2032: }
2033: }
2034: return 0;
2035: }
2037: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2038: {
2039: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2040: PetscInt dim,m,n,nrhs,minmn,maxmn;
2043: PetscFVGetSpatialDimension(fvm, &dim);
2044: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2045: ls->maxFaces = maxFaces;
2046: m = ls->maxFaces;
2047: n = dim;
2048: nrhs = ls->maxFaces;
2049: minmn = PetscMin(m,n);
2050: maxmn = PetscMax(m,n);
2051: ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2052: PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);
2053: return 0;
2054: }
2056: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2057: {
2058: fvm->ops->setfromoptions = NULL;
2059: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2060: fvm->ops->view = PetscFVView_LeastSquares;
2061: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2062: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2063: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2064: return 0;
2065: }
2067: /*MC
2068: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2070: Level: intermediate
2072: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2073: M*/
2075: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2076: {
2077: PetscFV_LeastSquares *ls;
2080: PetscNewLog(fvm, &ls);
2081: fvm->data = ls;
2083: ls->maxFaces = -1;
2084: ls->workSize = -1;
2085: ls->B = NULL;
2086: ls->Binv = NULL;
2087: ls->tau = NULL;
2088: ls->work = NULL;
2090: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2091: PetscFVInitialize_LeastSquares(fvm);
2092: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2093: return 0;
2094: }
2096: /*@
2097: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2099: Not collective
2101: Input parameters:
2102: + fvm - The PetscFV object
2103: - maxFaces - The maximum number of cell faces
2105: Level: intermediate
2107: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2108: @*/
2109: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2110: {
2112: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2113: return 0;
2114: }