Actual source code: fv.c
1: #include <petsc/private/petscfvimpl.h>
2: #include <petscdmplex.h>
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";
18: /*@C
19: PetscLimiterRegister - Adds a new PetscLimiter implementation
21: Not Collective
23: Input Parameters:
24: + name - The name of a new user-defined creation routine
25: - create_func - The creation routine itself
27: Notes:
28: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
30: Sample usage:
31: .vb
32: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
33: .ve
35: Then, your PetscLimiter type can be chosen with the procedural interface via
36: .vb
37: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
38: PetscLimiterSetType(PetscLimiter, "my_lim");
39: .ve
40: or at runtime via the option
41: .vb
42: -petsclimiter_type my_lim
43: .ve
45: Level: advanced
47: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
49: @*/
50: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
51: {
55: PetscFunctionListAdd(&PetscLimiterList, sname, function);
56: return(0);
57: }
59: /*@C
60: PetscLimiterSetType - Builds a particular PetscLimiter
62: Collective on lim
64: Input Parameters:
65: + lim - The PetscLimiter object
66: - name - The kind of limiter
68: Options Database Key:
69: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
71: Level: intermediate
73: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
74: @*/
75: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
76: {
77: PetscErrorCode (*r)(PetscLimiter);
78: PetscBool match;
83: PetscObjectTypeCompare((PetscObject) lim, name, &match);
84: if (match) return(0);
86: PetscLimiterRegisterAll();
87: PetscFunctionListFind(PetscLimiterList, name, &r);
88: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
90: if (lim->ops->destroy) {
91: (*lim->ops->destroy)(lim);
92: lim->ops->destroy = NULL;
93: }
94: (*r)(lim);
95: PetscObjectChangeTypeName((PetscObject) lim, name);
96: return(0);
97: }
99: /*@C
100: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
102: Not Collective
104: Input Parameter:
105: . lim - The PetscLimiter
107: Output Parameter:
108: . name - The PetscLimiter type name
110: Level: intermediate
112: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
113: @*/
114: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
115: {
121: PetscLimiterRegisterAll();
122: *name = ((PetscObject) lim)->type_name;
123: return(0);
124: }
126: /*@C
127: PetscLimiterViewFromOptions - View from Options
129: Collective on PetscLimiter
131: Input Parameters:
132: + A - the PetscLimiter object to view
133: . obj - Optional object
134: - name - command line option
136: Level: intermediate
137: .seealso: PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
138: @*/
139: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
140: {
145: PetscObjectViewFromOptions((PetscObject)A,obj,name);
146: return(0);
147: }
149: /*@C
150: PetscLimiterView - Views a PetscLimiter
152: Collective on lim
154: Input Parameter:
155: + lim - the PetscLimiter object to view
156: - v - the viewer
158: Level: beginner
160: .seealso: PetscLimiterDestroy()
161: @*/
162: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
163: {
168: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
169: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
170: return(0);
171: }
173: /*@
174: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
176: Collective on lim
178: Input Parameter:
179: . lim - the PetscLimiter object to set options for
181: Level: intermediate
183: .seealso: PetscLimiterView()
184: @*/
185: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
186: {
187: const char *defaultType;
188: char name[256];
189: PetscBool flg;
194: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
195: else defaultType = ((PetscObject) lim)->type_name;
196: PetscLimiterRegisterAll();
198: PetscObjectOptionsBegin((PetscObject) lim);
199: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
200: if (flg) {
201: PetscLimiterSetType(lim, name);
202: } else if (!((PetscObject) lim)->type_name) {
203: PetscLimiterSetType(lim, defaultType);
204: }
205: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
206: /* process any options handlers added with PetscObjectAddOptionsHandler() */
207: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
208: PetscOptionsEnd();
209: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
210: return(0);
211: }
213: /*@C
214: PetscLimiterSetUp - Construct data structures for the PetscLimiter
216: Collective on lim
218: Input Parameter:
219: . lim - the PetscLimiter object to setup
221: Level: intermediate
223: .seealso: PetscLimiterView(), PetscLimiterDestroy()
224: @*/
225: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
226: {
231: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
232: return(0);
233: }
235: /*@
236: PetscLimiterDestroy - Destroys a PetscLimiter object
238: Collective on lim
240: Input Parameter:
241: . lim - the PetscLimiter object to destroy
243: Level: beginner
245: .seealso: PetscLimiterView()
246: @*/
247: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
248: {
252: if (!*lim) return(0);
255: if (--((PetscObject)(*lim))->refct > 0) {*lim = NULL; return(0);}
256: ((PetscObject) (*lim))->refct = 0;
258: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
259: PetscHeaderDestroy(lim);
260: return(0);
261: }
263: /*@
264: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
266: Collective
268: Input Parameter:
269: . comm - The communicator for the PetscLimiter object
271: Output Parameter:
272: . lim - The PetscLimiter object
274: Level: beginner
276: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
277: @*/
278: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
279: {
280: PetscLimiter l;
285: PetscCitationsRegister(LimiterCitation,&Limitercite);
286: *lim = NULL;
287: PetscFVInitializePackage();
289: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
291: *lim = l;
292: return(0);
293: }
295: /*@
296: PetscLimiterLimit - Limit the flux
298: Input Parameters:
299: + lim - The PetscLimiter
300: - flim - The input field
302: Output Parameter:
303: . phi - The limited field
305: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
306: $ The classical flux-limited formulation is psi(r) where
307: $
308: $ r = (u[0] - u[-1]) / (u[1] - u[0])
309: $
310: $ The second order TVD region is bounded by
311: $
312: $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
313: $
314: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
315: $ phi(r)(r+1)/2 in which the reconstructed interface values are
316: $
317: $ u(v) = u[0] + phi(r) (grad u)[0] v
318: $
319: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
320: $
321: $ 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))
322: $
323: $ For a nicer symmetric formulation, rewrite in terms of
324: $
325: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
326: $
327: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
328: $
329: $ phi(r) = phi(1/r)
330: $
331: $ becomes
332: $
333: $ w(f) = w(1-f).
334: $
335: $ The limiters below implement this final form w(f). The reference methods are
336: $
337: $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
339: Level: beginner
341: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
342: @*/
343: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
344: {
350: (*lim->ops->limit)(lim, flim, phi);
351: return(0);
352: }
354: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
355: {
356: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
357: PetscErrorCode ierr;
360: PetscFree(l);
361: return(0);
362: }
364: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
365: {
366: PetscViewerFormat format;
367: PetscErrorCode ierr;
370: PetscViewerGetFormat(viewer, &format);
371: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
372: return(0);
373: }
375: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
376: {
377: PetscBool iascii;
383: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
384: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
385: return(0);
386: }
388: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
389: {
391: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
392: return(0);
393: }
395: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
396: {
398: lim->ops->view = PetscLimiterView_Sin;
399: lim->ops->destroy = PetscLimiterDestroy_Sin;
400: lim->ops->limit = PetscLimiterLimit_Sin;
401: return(0);
402: }
404: /*MC
405: PETSCLIMITERSIN = "sin" - A PetscLimiter object
407: Level: intermediate
409: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
410: M*/
412: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
413: {
414: PetscLimiter_Sin *l;
415: PetscErrorCode ierr;
419: PetscNewLog(lim, &l);
420: lim->data = l;
422: PetscLimiterInitialize_Sin(lim);
423: return(0);
424: }
426: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
427: {
428: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
429: PetscErrorCode ierr;
432: PetscFree(l);
433: return(0);
434: }
436: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
437: {
438: PetscViewerFormat format;
439: PetscErrorCode ierr;
442: PetscViewerGetFormat(viewer, &format);
443: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
444: return(0);
445: }
447: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
448: {
449: PetscBool iascii;
455: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
456: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
457: return(0);
458: }
460: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
461: {
463: *phi = 0.0;
464: return(0);
465: }
467: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
468: {
470: lim->ops->view = PetscLimiterView_Zero;
471: lim->ops->destroy = PetscLimiterDestroy_Zero;
472: lim->ops->limit = PetscLimiterLimit_Zero;
473: return(0);
474: }
476: /*MC
477: PETSCLIMITERZERO = "zero" - A PetscLimiter object
479: Level: intermediate
481: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
482: M*/
484: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
485: {
486: PetscLimiter_Zero *l;
487: PetscErrorCode ierr;
491: PetscNewLog(lim, &l);
492: lim->data = l;
494: PetscLimiterInitialize_Zero(lim);
495: return(0);
496: }
498: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
499: {
500: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
501: PetscErrorCode ierr;
504: PetscFree(l);
505: return(0);
506: }
508: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
509: {
510: PetscViewerFormat format;
511: PetscErrorCode ierr;
514: PetscViewerGetFormat(viewer, &format);
515: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
516: return(0);
517: }
519: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
520: {
521: PetscBool iascii;
527: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
528: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
529: return(0);
530: }
532: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
533: {
535: *phi = 1.0;
536: return(0);
537: }
539: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
540: {
542: lim->ops->view = PetscLimiterView_None;
543: lim->ops->destroy = PetscLimiterDestroy_None;
544: lim->ops->limit = PetscLimiterLimit_None;
545: return(0);
546: }
548: /*MC
549: PETSCLIMITERNONE = "none" - A PetscLimiter object
551: Level: intermediate
553: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
554: M*/
556: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
557: {
558: PetscLimiter_None *l;
559: PetscErrorCode ierr;
563: PetscNewLog(lim, &l);
564: lim->data = l;
566: PetscLimiterInitialize_None(lim);
567: return(0);
568: }
570: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
571: {
572: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
573: PetscErrorCode ierr;
576: PetscFree(l);
577: return(0);
578: }
580: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
581: {
582: PetscViewerFormat format;
583: PetscErrorCode ierr;
586: PetscViewerGetFormat(viewer, &format);
587: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
588: return(0);
589: }
591: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
592: {
593: PetscBool iascii;
599: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
600: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
601: return(0);
602: }
604: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
605: {
607: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
608: return(0);
609: }
611: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
612: {
614: lim->ops->view = PetscLimiterView_Minmod;
615: lim->ops->destroy = PetscLimiterDestroy_Minmod;
616: lim->ops->limit = PetscLimiterLimit_Minmod;
617: return(0);
618: }
620: /*MC
621: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
623: Level: intermediate
625: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
626: M*/
628: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
629: {
630: PetscLimiter_Minmod *l;
631: PetscErrorCode ierr;
635: PetscNewLog(lim, &l);
636: lim->data = l;
638: PetscLimiterInitialize_Minmod(lim);
639: return(0);
640: }
642: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
643: {
644: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
645: PetscErrorCode ierr;
648: PetscFree(l);
649: return(0);
650: }
652: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
653: {
654: PetscViewerFormat format;
655: PetscErrorCode ierr;
658: PetscViewerGetFormat(viewer, &format);
659: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
660: return(0);
661: }
663: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
664: {
665: PetscBool iascii;
671: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
672: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
673: return(0);
674: }
676: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
677: {
679: *phi = PetscMax(0, 4*f*(1-f));
680: return(0);
681: }
683: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
684: {
686: lim->ops->view = PetscLimiterView_VanLeer;
687: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
688: lim->ops->limit = PetscLimiterLimit_VanLeer;
689: return(0);
690: }
692: /*MC
693: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
695: Level: intermediate
697: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
698: M*/
700: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
701: {
702: PetscLimiter_VanLeer *l;
703: PetscErrorCode ierr;
707: PetscNewLog(lim, &l);
708: lim->data = l;
710: PetscLimiterInitialize_VanLeer(lim);
711: return(0);
712: }
714: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
715: {
716: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
717: PetscErrorCode ierr;
720: PetscFree(l);
721: return(0);
722: }
724: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
725: {
726: PetscViewerFormat format;
727: PetscErrorCode ierr;
730: PetscViewerGetFormat(viewer, &format);
731: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
732: return(0);
733: }
735: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
736: {
737: PetscBool iascii;
743: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
744: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
745: return(0);
746: }
748: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
749: {
751: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
752: return(0);
753: }
755: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
756: {
758: lim->ops->view = PetscLimiterView_VanAlbada;
759: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
760: lim->ops->limit = PetscLimiterLimit_VanAlbada;
761: return(0);
762: }
764: /*MC
765: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
767: Level: intermediate
769: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
770: M*/
772: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
773: {
774: PetscLimiter_VanAlbada *l;
775: PetscErrorCode ierr;
779: PetscNewLog(lim, &l);
780: lim->data = l;
782: PetscLimiterInitialize_VanAlbada(lim);
783: return(0);
784: }
786: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
787: {
788: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
789: PetscErrorCode ierr;
792: PetscFree(l);
793: return(0);
794: }
796: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
797: {
798: PetscViewerFormat format;
799: PetscErrorCode ierr;
802: PetscViewerGetFormat(viewer, &format);
803: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
804: return(0);
805: }
807: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
808: {
809: PetscBool iascii;
815: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
816: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
817: return(0);
818: }
820: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
821: {
823: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
824: return(0);
825: }
827: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
828: {
830: lim->ops->view = PetscLimiterView_Superbee;
831: lim->ops->destroy = PetscLimiterDestroy_Superbee;
832: lim->ops->limit = PetscLimiterLimit_Superbee;
833: return(0);
834: }
836: /*MC
837: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
839: Level: intermediate
841: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
842: M*/
844: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
845: {
846: PetscLimiter_Superbee *l;
847: PetscErrorCode ierr;
851: PetscNewLog(lim, &l);
852: lim->data = l;
854: PetscLimiterInitialize_Superbee(lim);
855: return(0);
856: }
858: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
859: {
860: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
861: PetscErrorCode ierr;
864: PetscFree(l);
865: return(0);
866: }
868: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
869: {
870: PetscViewerFormat format;
871: PetscErrorCode ierr;
874: PetscViewerGetFormat(viewer, &format);
875: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
876: return(0);
877: }
879: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
880: {
881: PetscBool iascii;
887: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
888: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
889: return(0);
890: }
892: /* aka Barth-Jespersen */
893: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
894: {
896: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
897: return(0);
898: }
900: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
901: {
903: lim->ops->view = PetscLimiterView_MC;
904: lim->ops->destroy = PetscLimiterDestroy_MC;
905: lim->ops->limit = PetscLimiterLimit_MC;
906: return(0);
907: }
909: /*MC
910: PETSCLIMITERMC = "mc" - A PetscLimiter object
912: Level: intermediate
914: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
915: M*/
917: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
918: {
919: PetscLimiter_MC *l;
920: PetscErrorCode ierr;
924: PetscNewLog(lim, &l);
925: lim->data = l;
927: PetscLimiterInitialize_MC(lim);
928: return(0);
929: }
931: PetscClassId PETSCFV_CLASSID = 0;
933: PetscFunctionList PetscFVList = NULL;
934: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
936: /*@C
937: PetscFVRegister - Adds a new PetscFV implementation
939: Not Collective
941: Input Parameters:
942: + name - The name of a new user-defined creation routine
943: - create_func - The creation routine itself
945: Notes:
946: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
948: Sample usage:
949: .vb
950: PetscFVRegister("my_fv", MyPetscFVCreate);
951: .ve
953: Then, your PetscFV type can be chosen with the procedural interface via
954: .vb
955: PetscFVCreate(MPI_Comm, PetscFV *);
956: PetscFVSetType(PetscFV, "my_fv");
957: .ve
958: or at runtime via the option
959: .vb
960: -petscfv_type my_fv
961: .ve
963: Level: advanced
965: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
967: @*/
968: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
969: {
973: PetscFunctionListAdd(&PetscFVList, sname, function);
974: return(0);
975: }
977: /*@C
978: PetscFVSetType - Builds a particular PetscFV
980: Collective on fvm
982: Input Parameters:
983: + fvm - The PetscFV object
984: - name - The kind of FVM space
986: Options Database Key:
987: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
989: Level: intermediate
991: .seealso: PetscFVGetType(), PetscFVCreate()
992: @*/
993: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
994: {
995: PetscErrorCode (*r)(PetscFV);
996: PetscBool match;
1001: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1002: if (match) return(0);
1004: PetscFVRegisterAll();
1005: PetscFunctionListFind(PetscFVList, name, &r);
1006: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1008: if (fvm->ops->destroy) {
1009: (*fvm->ops->destroy)(fvm);
1010: fvm->ops->destroy = NULL;
1011: }
1012: (*r)(fvm);
1013: PetscObjectChangeTypeName((PetscObject) fvm, name);
1014: return(0);
1015: }
1017: /*@C
1018: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1020: Not Collective
1022: Input Parameter:
1023: . fvm - The PetscFV
1025: Output Parameter:
1026: . name - The PetscFV type name
1028: Level: intermediate
1030: .seealso: PetscFVSetType(), PetscFVCreate()
1031: @*/
1032: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1033: {
1039: PetscFVRegisterAll();
1040: *name = ((PetscObject) fvm)->type_name;
1041: return(0);
1042: }
1044: /*@C
1045: PetscFVViewFromOptions - View from Options
1047: Collective on PetscFV
1049: Input Parameters:
1050: + A - the PetscFV object
1051: . obj - Optional object
1052: - name - command line option
1054: Level: intermediate
1055: .seealso: PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
1056: @*/
1057: PetscErrorCode PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
1058: {
1063: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1064: return(0);
1065: }
1067: /*@C
1068: PetscFVView - Views a PetscFV
1070: Collective on fvm
1072: Input Parameter:
1073: + fvm - the PetscFV object to view
1074: - v - the viewer
1076: Level: beginner
1078: .seealso: PetscFVDestroy()
1079: @*/
1080: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1081: {
1086: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1087: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1088: return(0);
1089: }
1091: /*@
1092: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1094: Collective on fvm
1096: Input Parameter:
1097: . fvm - the PetscFV object to set options for
1099: Options Database Key:
1100: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1102: Level: intermediate
1104: .seealso: PetscFVView()
1105: @*/
1106: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1107: {
1108: const char *defaultType;
1109: char name[256];
1110: PetscBool flg;
1115: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1116: else defaultType = ((PetscObject) fvm)->type_name;
1117: PetscFVRegisterAll();
1119: PetscObjectOptionsBegin((PetscObject) fvm);
1120: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1121: if (flg) {
1122: PetscFVSetType(fvm, name);
1123: } else if (!((PetscObject) fvm)->type_name) {
1124: PetscFVSetType(fvm, defaultType);
1126: }
1127: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1128: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1129: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1130: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1131: PetscLimiterSetFromOptions(fvm->limiter);
1132: PetscOptionsEnd();
1133: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1134: return(0);
1135: }
1137: /*@
1138: PetscFVSetUp - Construct data structures for the PetscFV
1140: Collective on fvm
1142: Input Parameter:
1143: . fvm - the PetscFV object to setup
1145: Level: intermediate
1147: .seealso: PetscFVView(), PetscFVDestroy()
1148: @*/
1149: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1150: {
1155: PetscLimiterSetUp(fvm->limiter);
1156: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1157: return(0);
1158: }
1160: /*@
1161: PetscFVDestroy - Destroys a PetscFV object
1163: Collective on fvm
1165: Input Parameter:
1166: . fvm - the PetscFV object to destroy
1168: Level: beginner
1170: .seealso: PetscFVView()
1171: @*/
1172: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1173: {
1174: PetscInt i;
1178: if (!*fvm) return(0);
1181: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; return(0);}
1182: ((PetscObject) (*fvm))->refct = 0;
1184: for (i = 0; i < (*fvm)->numComponents; i++) {
1185: PetscFree((*fvm)->componentNames[i]);
1186: }
1187: PetscFree((*fvm)->componentNames);
1188: PetscLimiterDestroy(&(*fvm)->limiter);
1189: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1190: PetscFree((*fvm)->fluxWork);
1191: PetscQuadratureDestroy(&(*fvm)->quadrature);
1192: PetscTabulationDestroy(&(*fvm)->T);
1194: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1195: PetscHeaderDestroy(fvm);
1196: return(0);
1197: }
1199: /*@
1200: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1202: Collective
1204: Input Parameter:
1205: . comm - The communicator for the PetscFV object
1207: Output Parameter:
1208: . fvm - The PetscFV object
1210: Level: beginner
1212: .seealso: PetscFVSetType(), PETSCFVUPWIND
1213: @*/
1214: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1215: {
1216: PetscFV f;
1221: *fvm = NULL;
1222: PetscFVInitializePackage();
1224: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1225: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1227: PetscLimiterCreate(comm, &f->limiter);
1228: f->numComponents = 1;
1229: f->dim = 0;
1230: f->computeGradients = PETSC_FALSE;
1231: f->fluxWork = NULL;
1232: PetscCalloc1(f->numComponents,&f->componentNames);
1234: *fvm = f;
1235: return(0);
1236: }
1238: /*@
1239: PetscFVSetLimiter - Set the limiter object
1241: Logically collective on fvm
1243: Input Parameters:
1244: + fvm - the PetscFV object
1245: - lim - The PetscLimiter
1247: Level: intermediate
1249: .seealso: PetscFVGetLimiter()
1250: @*/
1251: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1252: {
1258: PetscLimiterDestroy(&fvm->limiter);
1259: PetscObjectReference((PetscObject) lim);
1260: fvm->limiter = lim;
1261: return(0);
1262: }
1264: /*@
1265: PetscFVGetLimiter - Get the limiter object
1267: Not collective
1269: Input Parameter:
1270: . fvm - the PetscFV object
1272: Output Parameter:
1273: . lim - The PetscLimiter
1275: Level: intermediate
1277: .seealso: PetscFVSetLimiter()
1278: @*/
1279: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1280: {
1284: *lim = fvm->limiter;
1285: return(0);
1286: }
1288: /*@
1289: PetscFVSetNumComponents - Set the number of field components
1291: Logically collective on fvm
1293: Input Parameters:
1294: + fvm - the PetscFV object
1295: - comp - The number of components
1297: Level: intermediate
1299: .seealso: PetscFVGetNumComponents()
1300: @*/
1301: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1302: {
1307: if (fvm->numComponents != comp) {
1308: PetscInt i;
1310: for (i = 0; i < fvm->numComponents; i++) {
1311: PetscFree(fvm->componentNames[i]);
1312: }
1313: PetscFree(fvm->componentNames);
1314: PetscCalloc1(comp,&fvm->componentNames);
1315: }
1316: fvm->numComponents = comp;
1317: PetscFree(fvm->fluxWork);
1318: PetscMalloc1(comp, &fvm->fluxWork);
1319: return(0);
1320: }
1322: /*@
1323: PetscFVGetNumComponents - Get the number of field components
1325: Not collective
1327: Input Parameter:
1328: . fvm - the PetscFV object
1330: Output Parameter:
1331: , comp - The number of components
1333: Level: intermediate
1335: .seealso: PetscFVSetNumComponents()
1336: @*/
1337: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1338: {
1342: *comp = fvm->numComponents;
1343: return(0);
1344: }
1346: /*@C
1347: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1349: Logically collective on fvm
1350: Input Parameters:
1351: + fvm - the PetscFV object
1352: . comp - the component number
1353: - name - the component name
1355: Level: intermediate
1357: .seealso: PetscFVGetComponentName()
1358: @*/
1359: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1360: {
1364: PetscFree(fvm->componentNames[comp]);
1365: PetscStrallocpy(name,&fvm->componentNames[comp]);
1366: return(0);
1367: }
1369: /*@C
1370: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1372: Logically collective on fvm
1373: Input Parameters:
1374: + fvm - the PetscFV object
1375: - comp - the component number
1377: Output Parameter:
1378: . name - the component name
1380: Level: intermediate
1382: .seealso: PetscFVSetComponentName()
1383: @*/
1384: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1385: {
1387: *name = fvm->componentNames[comp];
1388: return(0);
1389: }
1391: /*@
1392: PetscFVSetSpatialDimension - Set the spatial dimension
1394: Logically collective on fvm
1396: Input Parameters:
1397: + fvm - the PetscFV object
1398: - dim - The spatial dimension
1400: Level: intermediate
1402: .seealso: PetscFVGetSpatialDimension()
1403: @*/
1404: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1405: {
1408: fvm->dim = dim;
1409: return(0);
1410: }
1412: /*@
1413: PetscFVGetSpatialDimension - Get the spatial dimension
1415: Logically collective on fvm
1417: Input Parameter:
1418: . fvm - the PetscFV object
1420: Output Parameter:
1421: . dim - The spatial dimension
1423: Level: intermediate
1425: .seealso: PetscFVSetSpatialDimension()
1426: @*/
1427: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1428: {
1432: *dim = fvm->dim;
1433: return(0);
1434: }
1436: /*@
1437: PetscFVSetComputeGradients - Toggle computation of cell gradients
1439: Logically collective on fvm
1441: Input Parameters:
1442: + fvm - the PetscFV object
1443: - computeGradients - Flag to compute cell gradients
1445: Level: intermediate
1447: .seealso: PetscFVGetComputeGradients()
1448: @*/
1449: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1450: {
1453: fvm->computeGradients = computeGradients;
1454: return(0);
1455: }
1457: /*@
1458: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1460: Not collective
1462: Input Parameter:
1463: . fvm - the PetscFV object
1465: Output Parameter:
1466: . computeGradients - Flag to compute cell gradients
1468: Level: intermediate
1470: .seealso: PetscFVSetComputeGradients()
1471: @*/
1472: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1473: {
1477: *computeGradients = fvm->computeGradients;
1478: return(0);
1479: }
1481: /*@
1482: PetscFVSetQuadrature - Set the quadrature object
1484: Logically collective on fvm
1486: Input Parameters:
1487: + fvm - the PetscFV object
1488: - q - The PetscQuadrature
1490: Level: intermediate
1492: .seealso: PetscFVGetQuadrature()
1493: @*/
1494: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1495: {
1500: PetscQuadratureDestroy(&fvm->quadrature);
1501: PetscObjectReference((PetscObject) q);
1502: fvm->quadrature = q;
1503: return(0);
1504: }
1506: /*@
1507: PetscFVGetQuadrature - Get the quadrature object
1509: Not collective
1511: Input Parameter:
1512: . fvm - the PetscFV object
1514: Output Parameter:
1515: . lim - The PetscQuadrature
1517: Level: intermediate
1519: .seealso: PetscFVSetQuadrature()
1520: @*/
1521: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1522: {
1526: if (!fvm->quadrature) {
1527: /* Create default 1-point quadrature */
1528: PetscReal *points, *weights;
1531: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1532: PetscCalloc1(fvm->dim, &points);
1533: PetscMalloc1(1, &weights);
1534: weights[0] = 1.0;
1535: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1536: }
1537: *q = fvm->quadrature;
1538: return(0);
1539: }
1541: /*@
1542: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1544: Not collective
1546: Input Parameter:
1547: . fvm - The PetscFV object
1549: Output Parameter:
1550: . sp - The PetscDualSpace object
1552: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1554: Level: intermediate
1556: .seealso: PetscFVCreate()
1557: @*/
1558: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1559: {
1563: if (!fvm->dualSpace) {
1564: DM K;
1565: PetscInt dim, Nc, c;
1566: PetscErrorCode ierr;
1568: PetscFVGetSpatialDimension(fvm, &dim);
1569: PetscFVGetNumComponents(fvm, &Nc);
1570: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1571: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1572: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1573: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1574: PetscDualSpaceSetDM(fvm->dualSpace, K);
1575: DMDestroy(&K);
1576: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1577: /* Should we be using PetscFVGetQuadrature() here? */
1578: for (c = 0; c < Nc; ++c) {
1579: PetscQuadrature qc;
1580: PetscReal *points, *weights;
1581: PetscErrorCode ierr;
1583: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1584: PetscCalloc1(dim, &points);
1585: PetscCalloc1(Nc, &weights);
1586: weights[c] = 1.0;
1587: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1588: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1589: PetscQuadratureDestroy(&qc);
1590: }
1591: PetscDualSpaceSetUp(fvm->dualSpace);
1592: }
1593: *sp = fvm->dualSpace;
1594: return(0);
1595: }
1597: /*@
1598: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1600: Not collective
1602: Input Parameters:
1603: + fvm - The PetscFV object
1604: - sp - The PetscDualSpace object
1606: Level: intermediate
1608: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1610: .seealso: PetscFVCreate()
1611: @*/
1612: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1613: {
1619: PetscDualSpaceDestroy(&fvm->dualSpace);
1620: fvm->dualSpace = sp;
1621: PetscObjectReference((PetscObject) fvm->dualSpace);
1622: return(0);
1623: }
1625: /*@C
1626: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1628: Not collective
1630: Input Parameter:
1631: . fvm - The PetscFV object
1633: Output Parameter:
1634: . T - The basis function values and derviatives at quadrature points
1636: Note:
1637: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1638: $ 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
1639: $ 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
1641: Level: intermediate
1643: .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1644: @*/
1645: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1646: {
1647: PetscInt npoints;
1648: const PetscReal *points;
1649: PetscErrorCode ierr;
1654: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1655: if (!fvm->T) {PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);}
1656: *T = fvm->T;
1657: return(0);
1658: }
1660: /*@C
1661: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1663: Not collective
1665: Input Parameters:
1666: + fvm - The PetscFV object
1667: . nrepl - The number of replicas
1668: . npoints - The number of tabulation points in a replica
1669: . points - The tabulation point coordinates
1670: - K - The order of derivative to tabulate
1672: Output Parameter:
1673: . T - The basis function values and derviative at tabulation points
1675: Note:
1676: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1677: $ 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
1678: $ 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
1680: Level: intermediate
1682: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1683: @*/
1684: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1685: {
1686: PetscInt pdim = 1; /* Dimension of approximation space P */
1687: PetscInt cdim; /* Spatial dimension */
1688: PetscInt Nc; /* Field components */
1689: PetscInt k, p, d, c, e;
1690: PetscErrorCode ierr;
1693: if (!npoints || K < 0) {
1694: *T = NULL;
1695: return(0);
1696: }
1700: PetscFVGetSpatialDimension(fvm, &cdim);
1701: PetscFVGetNumComponents(fvm, &Nc);
1702: PetscMalloc1(1, T);
1703: (*T)->K = !cdim ? 0 : K;
1704: (*T)->Nr = nrepl;
1705: (*T)->Np = npoints;
1706: (*T)->Nb = pdim;
1707: (*T)->Nc = Nc;
1708: (*T)->cdim = cdim;
1709: PetscMalloc1((*T)->K+1, &(*T)->T);
1710: for (k = 0; k <= (*T)->K; ++k) {
1711: PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
1712: }
1713: 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;}
1714: 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;}
1715: 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;}
1716: return(0);
1717: }
1719: /*@C
1720: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1722: Input Parameters:
1723: + fvm - The PetscFV object
1724: . numFaces - The number of cell faces which are not constrained
1725: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1727: Level: advanced
1729: .seealso: PetscFVCreate()
1730: @*/
1731: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1732: {
1737: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1738: return(0);
1739: }
1741: /*@C
1742: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1744: Not collective
1746: Input Parameters:
1747: + fvm - The PetscFV object for the field being integrated
1748: . prob - The PetscDS specifing the discretizations and continuum functions
1749: . field - The field being integrated
1750: . Nf - The number of faces in the chunk
1751: . fgeom - The face geometry for each face in the chunk
1752: . neighborVol - The volume for each pair of cells in the chunk
1753: . uL - The state from the cell on the left
1754: - uR - The state from the cell on the right
1756: Output Parameter:
1757: + fluxL - the left fluxes for each face
1758: - fluxR - the right fluxes for each face
1760: Level: developer
1762: .seealso: PetscFVCreate()
1763: @*/
1764: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1765: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1766: {
1771: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1772: return(0);
1773: }
1775: /*@
1776: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1777: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1778: sparsity). It is also used to create an interpolation between regularly refined meshes.
1780: Input Parameter:
1781: . fv - The initial PetscFV
1783: Output Parameter:
1784: . fvRef - The refined PetscFV
1786: Level: advanced
1788: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1789: @*/
1790: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1791: {
1792: PetscDualSpace Q, Qref;
1793: DM K, Kref;
1794: PetscQuadrature q, qref;
1795: DMPolytopeType ct;
1796: DMPlexCellRefiner cr;
1797: PetscReal *v0;
1798: PetscReal *jac, *invjac;
1799: PetscInt numComp, numSubelements, s;
1800: PetscErrorCode ierr;
1803: PetscFVGetDualSpace(fv, &Q);
1804: PetscFVGetQuadrature(fv, &q);
1805: PetscDualSpaceGetDM(Q, &K);
1806: /* Create dual space */
1807: PetscDualSpaceDuplicate(Q, &Qref);
1808: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1809: PetscDualSpaceSetDM(Qref, Kref);
1810: DMDestroy(&Kref);
1811: PetscDualSpaceSetUp(Qref);
1812: /* Create volume */
1813: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1814: PetscFVSetDualSpace(*fvRef, Qref);
1815: PetscFVGetNumComponents(fv, &numComp);
1816: PetscFVSetNumComponents(*fvRef, numComp);
1817: PetscFVSetUp(*fvRef);
1818: /* Create quadrature */
1819: DMPlexGetCellType(K, 0, &ct);
1820: DMPlexCellRefinerCreate(K, &cr);
1821: DMPlexCellRefinerGetAffineTransforms(cr, ct, &numSubelements, &v0, &jac, &invjac);
1822: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1823: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1824: for (s = 0; s < numSubelements; ++s) {
1825: PetscQuadrature qs;
1826: const PetscReal *points, *weights;
1827: PetscReal *p, *w;
1828: PetscInt dim, Nc, npoints, np;
1830: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1831: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1832: np = npoints/numSubelements;
1833: PetscMalloc1(np*dim,&p);
1834: PetscMalloc1(np*Nc,&w);
1835: PetscArraycpy(p, &points[s*np*dim], np*dim);
1836: PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1837: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1838: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1839: PetscQuadratureDestroy(&qs);
1840: }
1841: PetscFVSetQuadrature(*fvRef, qref);
1842: DMPlexCellRefinerDestroy(&cr);
1843: PetscQuadratureDestroy(&qref);
1844: PetscDualSpaceDestroy(&Qref);
1845: return(0);
1846: }
1848: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1849: {
1850: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1854: PetscFree(b);
1855: return(0);
1856: }
1858: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1859: {
1860: PetscInt Nc, c;
1861: PetscViewerFormat format;
1862: PetscErrorCode ierr;
1865: PetscFVGetNumComponents(fv, &Nc);
1866: PetscViewerGetFormat(viewer, &format);
1867: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1868: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1869: for (c = 0; c < Nc; c++) {
1870: if (fv->componentNames[c]) {
1871: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1872: }
1873: }
1874: return(0);
1875: }
1877: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1878: {
1879: PetscBool iascii;
1885: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1886: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1887: return(0);
1888: }
1890: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1891: {
1893: return(0);
1894: }
1896: /*
1897: neighborVol[f*2+0] contains the left geom
1898: neighborVol[f*2+1] contains the right geom
1899: */
1900: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1901: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1902: {
1903: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1904: void *rctx;
1905: PetscScalar *flux = fvm->fluxWork;
1906: const PetscScalar *constants;
1907: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1908: PetscErrorCode ierr;
1911: PetscDSGetTotalComponents(prob, &Nc);
1912: PetscDSGetTotalDimension(prob, &totDim);
1913: PetscDSGetFieldOffset(prob, field, &off);
1914: PetscDSGetRiemannSolver(prob, field, &riemann);
1915: PetscDSGetContext(prob, field, &rctx);
1916: PetscDSGetConstants(prob, &numConstants, &constants);
1917: PetscFVGetSpatialDimension(fvm, &dim);
1918: PetscFVGetNumComponents(fvm, &pdim);
1919: for (f = 0; f < Nf; ++f) {
1920: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1921: for (d = 0; d < pdim; ++d) {
1922: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1923: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1924: }
1925: }
1926: return(0);
1927: }
1929: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1930: {
1932: fvm->ops->setfromoptions = NULL;
1933: fvm->ops->setup = PetscFVSetUp_Upwind;
1934: fvm->ops->view = PetscFVView_Upwind;
1935: fvm->ops->destroy = PetscFVDestroy_Upwind;
1936: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1937: return(0);
1938: }
1940: /*MC
1941: PETSCFVUPWIND = "upwind" - A PetscFV object
1943: Level: intermediate
1945: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1946: M*/
1948: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1949: {
1950: PetscFV_Upwind *b;
1951: PetscErrorCode ierr;
1955: PetscNewLog(fvm,&b);
1956: fvm->data = b;
1958: PetscFVInitialize_Upwind(fvm);
1959: return(0);
1960: }
1962: #include <petscblaslapack.h>
1964: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1965: {
1966: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1967: PetscErrorCode ierr;
1970: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1971: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1972: PetscFree(ls);
1973: return(0);
1974: }
1976: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1977: {
1978: PetscInt Nc, c;
1979: PetscViewerFormat format;
1980: PetscErrorCode ierr;
1983: PetscFVGetNumComponents(fv, &Nc);
1984: PetscViewerGetFormat(viewer, &format);
1985: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1986: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1987: for (c = 0; c < Nc; c++) {
1988: if (fv->componentNames[c]) {
1989: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1990: }
1991: }
1992: return(0);
1993: }
1995: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1996: {
1997: PetscBool iascii;
2003: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2004: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2005: return(0);
2006: }
2008: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2009: {
2011: return(0);
2012: }
2014: /* Overwrites A. Can only handle full-rank problems with m>=n */
2015: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2016: {
2017: PetscBool debug = PETSC_FALSE;
2019: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2020: PetscScalar *R,*Q,*Aback,Alpha;
2023: if (debug) {
2024: PetscMalloc1(m*n,&Aback);
2025: PetscArraycpy(Aback,A,m*n);
2026: }
2028: PetscBLASIntCast(m,&M);
2029: PetscBLASIntCast(n,&N);
2030: PetscBLASIntCast(mstride,&lda);
2031: PetscBLASIntCast(worksize,&ldwork);
2032: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2033: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2034: PetscFPTrapPop();
2035: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2036: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2038: /* Extract an explicit representation of Q */
2039: Q = Ainv;
2040: PetscArraycpy(Q,A,mstride*n);
2041: K = N; /* full rank */
2042: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2043: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2045: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2046: Alpha = 1.0;
2047: ldb = lda;
2048: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2049: /* Ainv is Q, overwritten with inverse */
2051: if (debug) { /* Check that pseudo-inverse worked */
2052: PetscScalar Beta = 0.0;
2053: PetscBLASInt ldc;
2054: K = N;
2055: ldc = N;
2056: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2057: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2058: PetscFree(Aback);
2059: }
2060: return(0);
2061: }
2063: /* Overwrites A. Can handle degenerate problems and m<n. */
2064: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2065: {
2066: PetscBool debug = PETSC_FALSE;
2067: PetscScalar *Brhs, *Aback;
2068: PetscScalar *tmpwork;
2069: PetscReal rcond;
2070: #if defined (PETSC_USE_COMPLEX)
2071: PetscInt rworkSize;
2072: PetscReal *rwork;
2073: #endif
2074: PetscInt i, j, maxmn;
2075: PetscBLASInt M, N, lda, ldb, ldwork;
2076: PetscBLASInt nrhs, irank, info;
2080: if (debug) {
2081: PetscMalloc1(m*n,&Aback);
2082: PetscArraycpy(Aback,A,m*n);
2083: }
2085: /* initialize to identity */
2086: tmpwork = work;
2087: Brhs = Ainv;
2088: maxmn = PetscMax(m,n);
2089: for (j=0; j<maxmn; j++) {
2090: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2091: }
2093: PetscBLASIntCast(m,&M);
2094: PetscBLASIntCast(n,&N);
2095: PetscBLASIntCast(mstride,&lda);
2096: PetscBLASIntCast(maxmn,&ldb);
2097: PetscBLASIntCast(worksize,&ldwork);
2098: rcond = -1;
2099: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2100: nrhs = M;
2101: #if defined(PETSC_USE_COMPLEX)
2102: rworkSize = 5 * PetscMin(M,N);
2103: PetscMalloc1(rworkSize,&rwork);
2104: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info));
2105: PetscFPTrapPop();
2106: PetscFree(rwork);
2107: #else
2108: nrhs = M;
2109: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info));
2110: PetscFPTrapPop();
2111: #endif
2112: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2113: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2114: 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");
2115: return(0);
2116: }
2118: #if 0
2119: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2120: {
2121: PetscReal grad[2] = {0, 0};
2122: const PetscInt *faces;
2123: PetscInt numFaces, f;
2124: PetscErrorCode ierr;
2127: DMPlexGetConeSize(dm, cell, &numFaces);
2128: DMPlexGetCone(dm, cell, &faces);
2129: for (f = 0; f < numFaces; ++f) {
2130: const PetscInt *fcells;
2131: const CellGeom *cg1;
2132: const FaceGeom *fg;
2134: DMPlexGetSupport(dm, faces[f], &fcells);
2135: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2136: for (i = 0; i < 2; ++i) {
2137: PetscScalar du;
2139: if (fcells[i] == c) continue;
2140: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2141: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2142: grad[0] += fg->grad[!i][0] * du;
2143: grad[1] += fg->grad[!i][1] * du;
2144: }
2145: }
2146: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2147: return(0);
2148: }
2149: #endif
2151: /*
2152: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2154: Input Parameters:
2155: + fvm - The PetscFV object
2156: . numFaces - The number of cell faces which are not constrained
2157: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2159: Level: developer
2161: .seealso: PetscFVCreate()
2162: */
2163: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2164: {
2165: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2166: const PetscBool useSVD = PETSC_TRUE;
2167: const PetscInt maxFaces = ls->maxFaces;
2168: PetscInt dim, f, d;
2169: PetscErrorCode ierr;
2172: if (numFaces > maxFaces) {
2173: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2174: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2175: }
2176: PetscFVGetSpatialDimension(fvm, &dim);
2177: for (f = 0; f < numFaces; ++f) {
2178: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2179: }
2180: /* Overwrites B with garbage, returns Binv in row-major format */
2181: if (useSVD) {
2182: PetscInt maxmn = PetscMax(numFaces, dim);
2183: PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2184: /* Binv shaped in column-major, coldim=maxmn.*/
2185: for (f = 0; f < numFaces; ++f) {
2186: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d + maxmn*f];
2187: }
2188: } else {
2189: PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2190: /* Binv shaped in row-major, rowdim=maxFaces.*/
2191: for (f = 0; f < numFaces; ++f) {
2192: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces + f];
2193: }
2194: }
2195: return(0);
2196: }
2198: /*
2199: neighborVol[f*2+0] contains the left geom
2200: neighborVol[f*2+1] contains the right geom
2201: */
2202: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2203: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2204: {
2205: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2206: void *rctx;
2207: PetscScalar *flux = fvm->fluxWork;
2208: const PetscScalar *constants;
2209: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2210: PetscErrorCode ierr;
2213: PetscDSGetTotalComponents(prob, &Nc);
2214: PetscDSGetTotalDimension(prob, &totDim);
2215: PetscDSGetFieldOffset(prob, field, &off);
2216: PetscDSGetRiemannSolver(prob, field, &riemann);
2217: PetscDSGetContext(prob, field, &rctx);
2218: PetscDSGetConstants(prob, &numConstants, &constants);
2219: PetscFVGetSpatialDimension(fvm, &dim);
2220: PetscFVGetNumComponents(fvm, &pdim);
2221: for (f = 0; f < Nf; ++f) {
2222: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2223: for (d = 0; d < pdim; ++d) {
2224: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2225: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2226: }
2227: }
2228: return(0);
2229: }
2231: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2232: {
2233: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2234: PetscInt dim,m,n,nrhs,minmn,maxmn;
2235: PetscErrorCode ierr;
2239: PetscFVGetSpatialDimension(fvm, &dim);
2240: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2241: ls->maxFaces = maxFaces;
2242: m = ls->maxFaces;
2243: n = dim;
2244: nrhs = ls->maxFaces;
2245: minmn = PetscMin(m,n);
2246: maxmn = PetscMax(m,n);
2247: ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2248: PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);
2249: return(0);
2250: }
2252: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2253: {
2255: fvm->ops->setfromoptions = NULL;
2256: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2257: fvm->ops->view = PetscFVView_LeastSquares;
2258: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2259: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2260: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2261: return(0);
2262: }
2264: /*MC
2265: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2267: Level: intermediate
2269: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2270: M*/
2272: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2273: {
2274: PetscFV_LeastSquares *ls;
2275: PetscErrorCode ierr;
2279: PetscNewLog(fvm, &ls);
2280: fvm->data = ls;
2282: ls->maxFaces = -1;
2283: ls->workSize = -1;
2284: ls->B = NULL;
2285: ls->Binv = NULL;
2286: ls->tau = NULL;
2287: ls->work = NULL;
2289: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2290: PetscFVInitialize_LeastSquares(fvm);
2291: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2292: return(0);
2293: }
2295: /*@
2296: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2298: Not collective
2300: Input parameters:
2301: + fvm - The PetscFV object
2302: - maxFaces - The maximum number of cell faces
2304: Level: intermediate
2306: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2307: @*/
2308: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2309: {
2314: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2315: return(0);
2316: }