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: {
56: PetscFunctionListAdd(&PetscLimiterList, sname, function);
57: return(0);
58: }
60: /*@C
61: PetscLimiterSetType - Builds a particular PetscLimiter
63: Collective on lim
65: Input Parameters:
66: + lim - The PetscLimiter object
67: - name - The kind of limiter
69: Options Database Key:
70: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
72: Level: intermediate
74: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
75: @*/
76: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
77: {
78: PetscErrorCode (*r)(PetscLimiter);
79: PetscBool match;
84: PetscObjectTypeCompare((PetscObject) lim, name, &match);
85: if (match) return(0);
87: PetscLimiterRegisterAll();
88: PetscFunctionListFind(PetscLimiterList, name, &r);
89: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
91: if (lim->ops->destroy) {
92: (*lim->ops->destroy)(lim);
93: lim->ops->destroy = NULL;
94: }
95: (*r)(lim);
96: PetscObjectChangeTypeName((PetscObject) lim, name);
97: return(0);
98: }
100: /*@C
101: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
103: Not Collective
105: Input Parameter:
106: . lim - The PetscLimiter
108: Output Parameter:
109: . name - The PetscLimiter type name
111: Level: intermediate
113: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
114: @*/
115: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
116: {
122: PetscLimiterRegisterAll();
123: *name = ((PetscObject) lim)->type_name;
124: return(0);
125: }
127: /*@C
128: PetscLimiterViewFromOptions - View from Options
130: Collective on PetscLimiter
132: Input Parameters:
133: + A - the PetscLimiter object to view
134: . obj - Optional object
135: - name - command line option
137: Level: intermediate
138: .seealso: PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
139: @*/
140: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
141: {
146: PetscObjectViewFromOptions((PetscObject)A,obj,name);
147: return(0);
148: }
150: /*@C
151: PetscLimiterView - Views a PetscLimiter
153: Collective on lim
155: Input Parameters:
156: + lim - the PetscLimiter object to view
157: - v - the viewer
159: Level: beginner
161: .seealso: PetscLimiterDestroy()
162: @*/
163: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
164: {
169: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
170: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
171: return(0);
172: }
174: /*@
175: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
177: Collective on lim
179: Input Parameter:
180: . lim - the PetscLimiter object to set options for
182: Level: intermediate
184: .seealso: PetscLimiterView()
185: @*/
186: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
187: {
188: const char *defaultType;
189: char name[256];
190: PetscBool flg;
195: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
196: else defaultType = ((PetscObject) lim)->type_name;
197: PetscLimiterRegisterAll();
199: PetscObjectOptionsBegin((PetscObject) lim);
200: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
201: if (flg) {
202: PetscLimiterSetType(lim, name);
203: } else if (!((PetscObject) lim)->type_name) {
204: PetscLimiterSetType(lim, defaultType);
205: }
206: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
207: /* process any options handlers added with PetscObjectAddOptionsHandler() */
208: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
209: PetscOptionsEnd();
210: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
211: return(0);
212: }
214: /*@C
215: PetscLimiterSetUp - Construct data structures for the PetscLimiter
217: Collective on lim
219: Input Parameter:
220: . lim - the PetscLimiter object to setup
222: Level: intermediate
224: .seealso: PetscLimiterView(), PetscLimiterDestroy()
225: @*/
226: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
227: {
232: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
233: return(0);
234: }
236: /*@
237: PetscLimiterDestroy - Destroys a PetscLimiter object
239: Collective on lim
241: Input Parameter:
242: . lim - the PetscLimiter object to destroy
244: Level: beginner
246: .seealso: PetscLimiterView()
247: @*/
248: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
249: {
253: if (!*lim) return(0);
256: if (--((PetscObject)(*lim))->refct > 0) {*lim = NULL; return(0);}
257: ((PetscObject) (*lim))->refct = 0;
259: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
260: PetscHeaderDestroy(lim);
261: return(0);
262: }
264: /*@
265: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
267: Collective
269: Input Parameter:
270: . comm - The communicator for the PetscLimiter object
272: Output Parameter:
273: . lim - The PetscLimiter object
275: Level: beginner
277: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
278: @*/
279: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
280: {
281: PetscLimiter l;
286: PetscCitationsRegister(LimiterCitation,&Limitercite);
287: *lim = NULL;
288: PetscFVInitializePackage();
290: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
292: *lim = l;
293: return(0);
294: }
296: /*@
297: PetscLimiterLimit - Limit the flux
299: Input Parameters:
300: + lim - The PetscLimiter
301: - flim - The input field
303: Output Parameter:
304: . phi - The limited field
306: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
307: $ The classical flux-limited formulation is psi(r) where
308: $
309: $ r = (u[0] - u[-1]) / (u[1] - u[0])
310: $
311: $ The second order TVD region is bounded by
312: $
313: $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
314: $
315: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
316: $ phi(r)(r+1)/2 in which the reconstructed interface values are
317: $
318: $ u(v) = u[0] + phi(r) (grad u)[0] v
319: $
320: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
321: $
322: $ 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))
323: $
324: $ For a nicer symmetric formulation, rewrite in terms of
325: $
326: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
327: $
328: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
329: $
330: $ phi(r) = phi(1/r)
331: $
332: $ becomes
333: $
334: $ w(f) = w(1-f).
335: $
336: $ The limiters below implement this final form w(f). The reference methods are
337: $
338: $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
340: Level: beginner
342: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
343: @*/
344: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
345: {
351: (*lim->ops->limit)(lim, flim, phi);
352: return(0);
353: }
355: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
356: {
357: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
358: PetscErrorCode ierr;
361: PetscFree(l);
362: return(0);
363: }
365: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
366: {
367: PetscViewerFormat format;
368: PetscErrorCode ierr;
371: PetscViewerGetFormat(viewer, &format);
372: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
373: return(0);
374: }
376: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
377: {
378: PetscBool iascii;
384: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
385: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
386: return(0);
387: }
389: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
390: {
392: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
393: return(0);
394: }
396: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
397: {
399: lim->ops->view = PetscLimiterView_Sin;
400: lim->ops->destroy = PetscLimiterDestroy_Sin;
401: lim->ops->limit = PetscLimiterLimit_Sin;
402: return(0);
403: }
405: /*MC
406: PETSCLIMITERSIN = "sin" - A PetscLimiter object
408: Level: intermediate
410: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
411: M*/
413: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
414: {
415: PetscLimiter_Sin *l;
416: PetscErrorCode ierr;
420: PetscNewLog(lim, &l);
421: lim->data = l;
423: PetscLimiterInitialize_Sin(lim);
424: return(0);
425: }
427: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
428: {
429: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
430: PetscErrorCode ierr;
433: PetscFree(l);
434: return(0);
435: }
437: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
438: {
439: PetscViewerFormat format;
440: PetscErrorCode ierr;
443: PetscViewerGetFormat(viewer, &format);
444: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
445: return(0);
446: }
448: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
449: {
450: PetscBool iascii;
456: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
457: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
458: return(0);
459: }
461: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
462: {
464: *phi = 0.0;
465: return(0);
466: }
468: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
469: {
471: lim->ops->view = PetscLimiterView_Zero;
472: lim->ops->destroy = PetscLimiterDestroy_Zero;
473: lim->ops->limit = PetscLimiterLimit_Zero;
474: return(0);
475: }
477: /*MC
478: PETSCLIMITERZERO = "zero" - A PetscLimiter object
480: Level: intermediate
482: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
483: M*/
485: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
486: {
487: PetscLimiter_Zero *l;
488: PetscErrorCode ierr;
492: PetscNewLog(lim, &l);
493: lim->data = l;
495: PetscLimiterInitialize_Zero(lim);
496: return(0);
497: }
499: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
500: {
501: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
502: PetscErrorCode ierr;
505: PetscFree(l);
506: return(0);
507: }
509: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
510: {
511: PetscViewerFormat format;
512: PetscErrorCode ierr;
515: PetscViewerGetFormat(viewer, &format);
516: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
517: return(0);
518: }
520: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
521: {
522: PetscBool iascii;
528: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
529: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
530: return(0);
531: }
533: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
534: {
536: *phi = 1.0;
537: return(0);
538: }
540: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
541: {
543: lim->ops->view = PetscLimiterView_None;
544: lim->ops->destroy = PetscLimiterDestroy_None;
545: lim->ops->limit = PetscLimiterLimit_None;
546: return(0);
547: }
549: /*MC
550: PETSCLIMITERNONE = "none" - A PetscLimiter object
552: Level: intermediate
554: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
555: M*/
557: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
558: {
559: PetscLimiter_None *l;
560: PetscErrorCode ierr;
564: PetscNewLog(lim, &l);
565: lim->data = l;
567: PetscLimiterInitialize_None(lim);
568: return(0);
569: }
571: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
572: {
573: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
574: PetscErrorCode ierr;
577: PetscFree(l);
578: return(0);
579: }
581: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
582: {
583: PetscViewerFormat format;
584: PetscErrorCode ierr;
587: PetscViewerGetFormat(viewer, &format);
588: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
589: return(0);
590: }
592: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
593: {
594: PetscBool iascii;
600: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
601: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
602: return(0);
603: }
605: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
606: {
608: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
609: return(0);
610: }
612: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
613: {
615: lim->ops->view = PetscLimiterView_Minmod;
616: lim->ops->destroy = PetscLimiterDestroy_Minmod;
617: lim->ops->limit = PetscLimiterLimit_Minmod;
618: return(0);
619: }
621: /*MC
622: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
624: Level: intermediate
626: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
627: M*/
629: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
630: {
631: PetscLimiter_Minmod *l;
632: PetscErrorCode ierr;
636: PetscNewLog(lim, &l);
637: lim->data = l;
639: PetscLimiterInitialize_Minmod(lim);
640: return(0);
641: }
643: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
644: {
645: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
646: PetscErrorCode ierr;
649: PetscFree(l);
650: return(0);
651: }
653: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
654: {
655: PetscViewerFormat format;
656: PetscErrorCode ierr;
659: PetscViewerGetFormat(viewer, &format);
660: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
661: return(0);
662: }
664: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
665: {
666: PetscBool iascii;
672: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
673: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
674: return(0);
675: }
677: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
678: {
680: *phi = PetscMax(0, 4*f*(1-f));
681: return(0);
682: }
684: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
685: {
687: lim->ops->view = PetscLimiterView_VanLeer;
688: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
689: lim->ops->limit = PetscLimiterLimit_VanLeer;
690: return(0);
691: }
693: /*MC
694: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
696: Level: intermediate
698: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
699: M*/
701: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
702: {
703: PetscLimiter_VanLeer *l;
704: PetscErrorCode ierr;
708: PetscNewLog(lim, &l);
709: lim->data = l;
711: PetscLimiterInitialize_VanLeer(lim);
712: return(0);
713: }
715: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
716: {
717: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
718: PetscErrorCode ierr;
721: PetscFree(l);
722: return(0);
723: }
725: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
726: {
727: PetscViewerFormat format;
728: PetscErrorCode ierr;
731: PetscViewerGetFormat(viewer, &format);
732: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
733: return(0);
734: }
736: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
737: {
738: PetscBool iascii;
744: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
745: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
746: return(0);
747: }
749: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
750: {
752: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
753: return(0);
754: }
756: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
757: {
759: lim->ops->view = PetscLimiterView_VanAlbada;
760: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
761: lim->ops->limit = PetscLimiterLimit_VanAlbada;
762: return(0);
763: }
765: /*MC
766: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
768: Level: intermediate
770: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
771: M*/
773: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
774: {
775: PetscLimiter_VanAlbada *l;
776: PetscErrorCode ierr;
780: PetscNewLog(lim, &l);
781: lim->data = l;
783: PetscLimiterInitialize_VanAlbada(lim);
784: return(0);
785: }
787: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
788: {
789: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
790: PetscErrorCode ierr;
793: PetscFree(l);
794: return(0);
795: }
797: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
798: {
799: PetscViewerFormat format;
800: PetscErrorCode ierr;
803: PetscViewerGetFormat(viewer, &format);
804: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
805: return(0);
806: }
808: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
809: {
810: PetscBool iascii;
816: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
817: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
818: return(0);
819: }
821: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
822: {
824: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
825: return(0);
826: }
828: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
829: {
831: lim->ops->view = PetscLimiterView_Superbee;
832: lim->ops->destroy = PetscLimiterDestroy_Superbee;
833: lim->ops->limit = PetscLimiterLimit_Superbee;
834: return(0);
835: }
837: /*MC
838: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
840: Level: intermediate
842: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
843: M*/
845: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
846: {
847: PetscLimiter_Superbee *l;
848: PetscErrorCode ierr;
852: PetscNewLog(lim, &l);
853: lim->data = l;
855: PetscLimiterInitialize_Superbee(lim);
856: return(0);
857: }
859: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
860: {
861: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
862: PetscErrorCode ierr;
865: PetscFree(l);
866: return(0);
867: }
869: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
870: {
871: PetscViewerFormat format;
872: PetscErrorCode ierr;
875: PetscViewerGetFormat(viewer, &format);
876: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
877: return(0);
878: }
880: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
881: {
882: PetscBool iascii;
888: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
889: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
890: return(0);
891: }
893: /* aka Barth-Jespersen */
894: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
895: {
897: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
898: return(0);
899: }
901: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
902: {
904: lim->ops->view = PetscLimiterView_MC;
905: lim->ops->destroy = PetscLimiterDestroy_MC;
906: lim->ops->limit = PetscLimiterLimit_MC;
907: return(0);
908: }
910: /*MC
911: PETSCLIMITERMC = "mc" - A PetscLimiter object
913: Level: intermediate
915: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
916: M*/
918: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
919: {
920: PetscLimiter_MC *l;
921: PetscErrorCode ierr;
925: PetscNewLog(lim, &l);
926: lim->data = l;
928: PetscLimiterInitialize_MC(lim);
929: return(0);
930: }
932: PetscClassId PETSCFV_CLASSID = 0;
934: PetscFunctionList PetscFVList = NULL;
935: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
937: /*@C
938: PetscFVRegister - Adds a new PetscFV implementation
940: Not Collective
942: Input Parameters:
943: + name - The name of a new user-defined creation routine
944: - create_func - The creation routine itself
946: Notes:
947: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
949: Sample usage:
950: .vb
951: PetscFVRegister("my_fv", MyPetscFVCreate);
952: .ve
954: Then, your PetscFV type can be chosen with the procedural interface via
955: .vb
956: PetscFVCreate(MPI_Comm, PetscFV *);
957: PetscFVSetType(PetscFV, "my_fv");
958: .ve
959: or at runtime via the option
960: .vb
961: -petscfv_type my_fv
962: .ve
964: Level: advanced
966: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
968: @*/
969: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
970: {
974: PetscFunctionListAdd(&PetscFVList, sname, function);
975: return(0);
976: }
978: /*@C
979: PetscFVSetType - Builds a particular PetscFV
981: Collective on fvm
983: Input Parameters:
984: + fvm - The PetscFV object
985: - name - The kind of FVM space
987: Options Database Key:
988: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
990: Level: intermediate
992: .seealso: PetscFVGetType(), PetscFVCreate()
993: @*/
994: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
995: {
996: PetscErrorCode (*r)(PetscFV);
997: PetscBool match;
1002: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1003: if (match) return(0);
1005: PetscFVRegisterAll();
1006: PetscFunctionListFind(PetscFVList, name, &r);
1007: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1009: if (fvm->ops->destroy) {
1010: (*fvm->ops->destroy)(fvm);
1011: fvm->ops->destroy = NULL;
1012: }
1013: (*r)(fvm);
1014: PetscObjectChangeTypeName((PetscObject) fvm, name);
1015: return(0);
1016: }
1018: /*@C
1019: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1021: Not Collective
1023: Input Parameter:
1024: . fvm - The PetscFV
1026: Output Parameter:
1027: . name - The PetscFV type name
1029: Level: intermediate
1031: .seealso: PetscFVSetType(), PetscFVCreate()
1032: @*/
1033: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1034: {
1040: PetscFVRegisterAll();
1041: *name = ((PetscObject) fvm)->type_name;
1042: return(0);
1043: }
1045: /*@C
1046: PetscFVViewFromOptions - View from Options
1048: Collective on PetscFV
1050: Input Parameters:
1051: + A - the PetscFV object
1052: . obj - Optional object
1053: - name - command line option
1055: Level: intermediate
1056: .seealso: PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
1057: @*/
1058: PetscErrorCode PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
1059: {
1064: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1065: return(0);
1066: }
1068: /*@C
1069: PetscFVView - Views a PetscFV
1071: Collective on fvm
1073: Input Parameters:
1074: + fvm - the PetscFV object to view
1075: - v - the viewer
1077: Level: beginner
1079: .seealso: PetscFVDestroy()
1080: @*/
1081: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1082: {
1087: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1088: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1089: return(0);
1090: }
1092: /*@
1093: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1095: Collective on fvm
1097: Input Parameter:
1098: . fvm - the PetscFV object to set options for
1100: Options Database Key:
1101: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1103: Level: intermediate
1105: .seealso: PetscFVView()
1106: @*/
1107: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1108: {
1109: const char *defaultType;
1110: char name[256];
1111: PetscBool flg;
1116: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1117: else defaultType = ((PetscObject) fvm)->type_name;
1118: PetscFVRegisterAll();
1120: PetscObjectOptionsBegin((PetscObject) fvm);
1121: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1122: if (flg) {
1123: PetscFVSetType(fvm, name);
1124: } else if (!((PetscObject) fvm)->type_name) {
1125: PetscFVSetType(fvm, defaultType);
1127: }
1128: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1129: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1130: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1131: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1132: PetscLimiterSetFromOptions(fvm->limiter);
1133: PetscOptionsEnd();
1134: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1135: return(0);
1136: }
1138: /*@
1139: PetscFVSetUp - Construct data structures for the PetscFV
1141: Collective on fvm
1143: Input Parameter:
1144: . fvm - the PetscFV object to setup
1146: Level: intermediate
1148: .seealso: PetscFVView(), PetscFVDestroy()
1149: @*/
1150: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1151: {
1156: PetscLimiterSetUp(fvm->limiter);
1157: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1158: return(0);
1159: }
1161: /*@
1162: PetscFVDestroy - Destroys a PetscFV object
1164: Collective on fvm
1166: Input Parameter:
1167: . fvm - the PetscFV object to destroy
1169: Level: beginner
1171: .seealso: PetscFVView()
1172: @*/
1173: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1174: {
1175: PetscInt i;
1179: if (!*fvm) return(0);
1182: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; return(0);}
1183: ((PetscObject) (*fvm))->refct = 0;
1185: for (i = 0; i < (*fvm)->numComponents; i++) {
1186: PetscFree((*fvm)->componentNames[i]);
1187: }
1188: PetscFree((*fvm)->componentNames);
1189: PetscLimiterDestroy(&(*fvm)->limiter);
1190: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1191: PetscFree((*fvm)->fluxWork);
1192: PetscQuadratureDestroy(&(*fvm)->quadrature);
1193: PetscTabulationDestroy(&(*fvm)->T);
1195: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1196: PetscHeaderDestroy(fvm);
1197: return(0);
1198: }
1200: /*@
1201: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1203: Collective
1205: Input Parameter:
1206: . comm - The communicator for the PetscFV object
1208: Output Parameter:
1209: . fvm - The PetscFV object
1211: Level: beginner
1213: .seealso: PetscFVSetType(), PETSCFVUPWIND
1214: @*/
1215: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1216: {
1217: PetscFV f;
1222: *fvm = NULL;
1223: PetscFVInitializePackage();
1225: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1226: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1228: PetscLimiterCreate(comm, &f->limiter);
1229: f->numComponents = 1;
1230: f->dim = 0;
1231: f->computeGradients = PETSC_FALSE;
1232: f->fluxWork = NULL;
1233: PetscCalloc1(f->numComponents,&f->componentNames);
1235: *fvm = f;
1236: return(0);
1237: }
1239: /*@
1240: PetscFVSetLimiter - Set the limiter object
1242: Logically collective on fvm
1244: Input Parameters:
1245: + fvm - the PetscFV object
1246: - lim - The PetscLimiter
1248: Level: intermediate
1250: .seealso: PetscFVGetLimiter()
1251: @*/
1252: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1253: {
1259: PetscLimiterDestroy(&fvm->limiter);
1260: PetscObjectReference((PetscObject) lim);
1261: fvm->limiter = lim;
1262: return(0);
1263: }
1265: /*@
1266: PetscFVGetLimiter - Get the limiter object
1268: Not collective
1270: Input Parameter:
1271: . fvm - the PetscFV object
1273: Output Parameter:
1274: . lim - The PetscLimiter
1276: Level: intermediate
1278: .seealso: PetscFVSetLimiter()
1279: @*/
1280: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1281: {
1285: *lim = fvm->limiter;
1286: return(0);
1287: }
1289: /*@
1290: PetscFVSetNumComponents - Set the number of field components
1292: Logically collective on fvm
1294: Input Parameters:
1295: + fvm - the PetscFV object
1296: - comp - The number of components
1298: Level: intermediate
1300: .seealso: PetscFVGetNumComponents()
1301: @*/
1302: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1303: {
1308: if (fvm->numComponents != comp) {
1309: PetscInt i;
1311: for (i = 0; i < fvm->numComponents; i++) {
1312: PetscFree(fvm->componentNames[i]);
1313: }
1314: PetscFree(fvm->componentNames);
1315: PetscCalloc1(comp,&fvm->componentNames);
1316: }
1317: fvm->numComponents = comp;
1318: PetscFree(fvm->fluxWork);
1319: PetscMalloc1(comp, &fvm->fluxWork);
1320: return(0);
1321: }
1323: /*@
1324: PetscFVGetNumComponents - Get the number of field components
1326: Not collective
1328: Input Parameter:
1329: . fvm - the PetscFV object
1331: Output Parameter:
1332: , comp - The number of components
1334: Level: intermediate
1336: .seealso: PetscFVSetNumComponents()
1337: @*/
1338: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1339: {
1343: *comp = fvm->numComponents;
1344: return(0);
1345: }
1347: /*@C
1348: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1350: Logically collective on fvm
1351: Input Parameters:
1352: + fvm - the PetscFV object
1353: . comp - the component number
1354: - name - the component name
1356: Level: intermediate
1358: .seealso: PetscFVGetComponentName()
1359: @*/
1360: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1361: {
1365: PetscFree(fvm->componentNames[comp]);
1366: PetscStrallocpy(name,&fvm->componentNames[comp]);
1367: return(0);
1368: }
1370: /*@C
1371: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1373: Logically collective on fvm
1374: Input Parameters:
1375: + fvm - the PetscFV object
1376: - comp - the component number
1378: Output Parameter:
1379: . name - the component name
1381: Level: intermediate
1383: .seealso: PetscFVSetComponentName()
1384: @*/
1385: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1386: {
1388: *name = fvm->componentNames[comp];
1389: return(0);
1390: }
1392: /*@
1393: PetscFVSetSpatialDimension - Set the spatial dimension
1395: Logically collective on fvm
1397: Input Parameters:
1398: + fvm - the PetscFV object
1399: - dim - The spatial dimension
1401: Level: intermediate
1403: .seealso: PetscFVGetSpatialDimension()
1404: @*/
1405: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1406: {
1409: fvm->dim = dim;
1410: return(0);
1411: }
1413: /*@
1414: PetscFVGetSpatialDimension - Get the spatial dimension
1416: Logically collective on fvm
1418: Input Parameter:
1419: . fvm - the PetscFV object
1421: Output Parameter:
1422: . dim - The spatial dimension
1424: Level: intermediate
1426: .seealso: PetscFVSetSpatialDimension()
1427: @*/
1428: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1429: {
1433: *dim = fvm->dim;
1434: return(0);
1435: }
1437: /*@
1438: PetscFVSetComputeGradients - Toggle computation of cell gradients
1440: Logically collective on fvm
1442: Input Parameters:
1443: + fvm - the PetscFV object
1444: - computeGradients - Flag to compute cell gradients
1446: Level: intermediate
1448: .seealso: PetscFVGetComputeGradients()
1449: @*/
1450: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1451: {
1454: fvm->computeGradients = computeGradients;
1455: return(0);
1456: }
1458: /*@
1459: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1461: Not collective
1463: Input Parameter:
1464: . fvm - the PetscFV object
1466: Output Parameter:
1467: . computeGradients - Flag to compute cell gradients
1469: Level: intermediate
1471: .seealso: PetscFVSetComputeGradients()
1472: @*/
1473: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1474: {
1478: *computeGradients = fvm->computeGradients;
1479: return(0);
1480: }
1482: /*@
1483: PetscFVSetQuadrature - Set the quadrature object
1485: Logically collective on fvm
1487: Input Parameters:
1488: + fvm - the PetscFV object
1489: - q - The PetscQuadrature
1491: Level: intermediate
1493: .seealso: PetscFVGetQuadrature()
1494: @*/
1495: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1496: {
1501: PetscQuadratureDestroy(&fvm->quadrature);
1502: PetscObjectReference((PetscObject) q);
1503: fvm->quadrature = q;
1504: return(0);
1505: }
1507: /*@
1508: PetscFVGetQuadrature - Get the quadrature object
1510: Not collective
1512: Input Parameter:
1513: . fvm - the PetscFV object
1515: Output Parameter:
1516: . lim - The PetscQuadrature
1518: Level: intermediate
1520: .seealso: PetscFVSetQuadrature()
1521: @*/
1522: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1523: {
1527: if (!fvm->quadrature) {
1528: /* Create default 1-point quadrature */
1529: PetscReal *points, *weights;
1532: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1533: PetscCalloc1(fvm->dim, &points);
1534: PetscMalloc1(1, &weights);
1535: weights[0] = 1.0;
1536: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1537: }
1538: *q = fvm->quadrature;
1539: return(0);
1540: }
1542: /*@
1543: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1545: Not collective
1547: Input Parameter:
1548: . fvm - The PetscFV object
1550: Output Parameter:
1551: . sp - The PetscDualSpace object
1553: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1555: Level: intermediate
1557: .seealso: PetscFVCreate()
1558: @*/
1559: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1560: {
1564: if (!fvm->dualSpace) {
1565: DM K;
1566: PetscInt dim, Nc, c;
1567: PetscErrorCode ierr;
1569: PetscFVGetSpatialDimension(fvm, &dim);
1570: PetscFVGetNumComponents(fvm, &Nc);
1571: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1572: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1573: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1574: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1575: PetscDualSpaceSetDM(fvm->dualSpace, K);
1576: DMDestroy(&K);
1577: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1578: /* Should we be using PetscFVGetQuadrature() here? */
1579: for (c = 0; c < Nc; ++c) {
1580: PetscQuadrature qc;
1581: PetscReal *points, *weights;
1582: PetscErrorCode ierr;
1584: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1585: PetscCalloc1(dim, &points);
1586: PetscCalloc1(Nc, &weights);
1587: weights[c] = 1.0;
1588: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1589: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1590: PetscQuadratureDestroy(&qc);
1591: }
1592: PetscDualSpaceSetUp(fvm->dualSpace);
1593: }
1594: *sp = fvm->dualSpace;
1595: return(0);
1596: }
1598: /*@
1599: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1601: Not collective
1603: Input Parameters:
1604: + fvm - The PetscFV object
1605: - sp - The PetscDualSpace object
1607: Level: intermediate
1609: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1611: .seealso: PetscFVCreate()
1612: @*/
1613: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1614: {
1620: PetscDualSpaceDestroy(&fvm->dualSpace);
1621: fvm->dualSpace = sp;
1622: PetscObjectReference((PetscObject) fvm->dualSpace);
1623: return(0);
1624: }
1626: /*@C
1627: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1629: Not collective
1631: Input Parameter:
1632: . fvm - The PetscFV object
1634: Output Parameter:
1635: . T - The basis function values and derivatives at quadrature points
1637: Note:
1638: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1639: $ 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
1640: $ 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
1642: Level: intermediate
1644: .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1645: @*/
1646: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1647: {
1648: PetscInt npoints;
1649: const PetscReal *points;
1650: PetscErrorCode ierr;
1655: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1656: if (!fvm->T) {PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);}
1657: *T = fvm->T;
1658: return(0);
1659: }
1661: /*@C
1662: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1664: Not collective
1666: Input Parameters:
1667: + fvm - The PetscFV object
1668: . nrepl - The number of replicas
1669: . npoints - The number of tabulation points in a replica
1670: . points - The tabulation point coordinates
1671: - K - The order of derivative to tabulate
1673: Output Parameter:
1674: . T - The basis function values and derivative at tabulation points
1676: Note:
1677: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1678: $ 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
1679: $ 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
1681: Level: intermediate
1683: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1684: @*/
1685: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1686: {
1687: PetscInt pdim = 1; /* Dimension of approximation space P */
1688: PetscInt cdim; /* Spatial dimension */
1689: PetscInt Nc; /* Field components */
1690: PetscInt k, p, d, c, e;
1691: PetscErrorCode ierr;
1694: if (!npoints || K < 0) {
1695: *T = NULL;
1696: return(0);
1697: }
1701: PetscFVGetSpatialDimension(fvm, &cdim);
1702: PetscFVGetNumComponents(fvm, &Nc);
1703: PetscMalloc1(1, T);
1704: (*T)->K = !cdim ? 0 : K;
1705: (*T)->Nr = nrepl;
1706: (*T)->Np = npoints;
1707: (*T)->Nb = pdim;
1708: (*T)->Nc = Nc;
1709: (*T)->cdim = cdim;
1710: PetscMalloc1((*T)->K+1, &(*T)->T);
1711: for (k = 0; k <= (*T)->K; ++k) {
1712: PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
1713: }
1714: 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;}
1715: 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;}
1716: 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;}
1717: return(0);
1718: }
1720: /*@C
1721: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1723: Input Parameters:
1724: + fvm - The PetscFV object
1725: . numFaces - The number of cell faces which are not constrained
1726: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1728: Level: advanced
1730: .seealso: PetscFVCreate()
1731: @*/
1732: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1733: {
1738: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1739: return(0);
1740: }
1742: /*@C
1743: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1745: Not collective
1747: Input Parameters:
1748: + fvm - The PetscFV object for the field being integrated
1749: . prob - The PetscDS specifing the discretizations and continuum functions
1750: . field - The field being integrated
1751: . Nf - The number of faces in the chunk
1752: . fgeom - The face geometry for each face in the chunk
1753: . neighborVol - The volume for each pair of cells in the chunk
1754: . uL - The state from the cell on the left
1755: - uR - The state from the cell on the right
1757: Output Parameters:
1758: + fluxL - the left fluxes for each face
1759: - fluxR - the right fluxes for each face
1761: Level: developer
1763: .seealso: PetscFVCreate()
1764: @*/
1765: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1766: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1767: {
1772: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1773: return(0);
1774: }
1776: /*@
1777: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1778: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1779: sparsity). It is also used to create an interpolation between regularly refined meshes.
1781: Input Parameter:
1782: . fv - The initial PetscFV
1784: Output Parameter:
1785: . fvRef - The refined PetscFV
1787: Level: advanced
1789: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1790: @*/
1791: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1792: {
1793: PetscDualSpace Q, Qref;
1794: DM K, Kref;
1795: PetscQuadrature q, qref;
1796: DMPolytopeType ct;
1797: DMPlexTransform tr;
1798: PetscReal *v0;
1799: PetscReal *jac, *invjac;
1800: PetscInt numComp, numSubelements, s;
1801: PetscErrorCode ierr;
1804: PetscFVGetDualSpace(fv, &Q);
1805: PetscFVGetQuadrature(fv, &q);
1806: PetscDualSpaceGetDM(Q, &K);
1807: /* Create dual space */
1808: PetscDualSpaceDuplicate(Q, &Qref);
1809: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1810: PetscDualSpaceSetDM(Qref, Kref);
1811: DMDestroy(&Kref);
1812: PetscDualSpaceSetUp(Qref);
1813: /* Create volume */
1814: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1815: PetscFVSetDualSpace(*fvRef, Qref);
1816: PetscFVGetNumComponents(fv, &numComp);
1817: PetscFVSetNumComponents(*fvRef, numComp);
1818: PetscFVSetUp(*fvRef);
1819: /* Create quadrature */
1820: DMPlexGetCellType(K, 0, &ct);
1821: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1822: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1823: DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1824: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1825: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1826: for (s = 0; s < numSubelements; ++s) {
1827: PetscQuadrature qs;
1828: const PetscReal *points, *weights;
1829: PetscReal *p, *w;
1830: PetscInt dim, Nc, npoints, np;
1832: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1833: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1834: np = npoints/numSubelements;
1835: PetscMalloc1(np*dim,&p);
1836: PetscMalloc1(np*Nc,&w);
1837: PetscArraycpy(p, &points[s*np*dim], np*dim);
1838: PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1839: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1840: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1841: PetscQuadratureDestroy(&qs);
1842: }
1843: PetscFVSetQuadrature(*fvRef, qref);
1844: DMPlexTransformDestroy(&tr);
1845: PetscQuadratureDestroy(&qref);
1846: PetscDualSpaceDestroy(&Qref);
1847: return(0);
1848: }
1850: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1851: {
1852: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1856: PetscFree(b);
1857: return(0);
1858: }
1860: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1861: {
1862: PetscInt Nc, c;
1863: PetscViewerFormat format;
1864: PetscErrorCode ierr;
1867: PetscFVGetNumComponents(fv, &Nc);
1868: PetscViewerGetFormat(viewer, &format);
1869: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1870: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1871: for (c = 0; c < Nc; c++) {
1872: if (fv->componentNames[c]) {
1873: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1874: }
1875: }
1876: return(0);
1877: }
1879: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1880: {
1881: PetscBool iascii;
1887: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1888: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1889: return(0);
1890: }
1892: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1893: {
1895: return(0);
1896: }
1898: /*
1899: neighborVol[f*2+0] contains the left geom
1900: neighborVol[f*2+1] contains the right geom
1901: */
1902: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1903: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1904: {
1905: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1906: void *rctx;
1907: PetscScalar *flux = fvm->fluxWork;
1908: const PetscScalar *constants;
1909: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1910: PetscErrorCode ierr;
1913: PetscDSGetTotalComponents(prob, &Nc);
1914: PetscDSGetTotalDimension(prob, &totDim);
1915: PetscDSGetFieldOffset(prob, field, &off);
1916: PetscDSGetRiemannSolver(prob, field, &riemann);
1917: PetscDSGetContext(prob, field, &rctx);
1918: PetscDSGetConstants(prob, &numConstants, &constants);
1919: PetscFVGetSpatialDimension(fvm, &dim);
1920: PetscFVGetNumComponents(fvm, &pdim);
1921: for (f = 0; f < Nf; ++f) {
1922: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1923: for (d = 0; d < pdim; ++d) {
1924: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1925: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1926: }
1927: }
1928: return(0);
1929: }
1931: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1932: {
1934: fvm->ops->setfromoptions = NULL;
1935: fvm->ops->setup = PetscFVSetUp_Upwind;
1936: fvm->ops->view = PetscFVView_Upwind;
1937: fvm->ops->destroy = PetscFVDestroy_Upwind;
1938: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1939: return(0);
1940: }
1942: /*MC
1943: PETSCFVUPWIND = "upwind" - A PetscFV object
1945: Level: intermediate
1947: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1948: M*/
1950: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1951: {
1952: PetscFV_Upwind *b;
1953: PetscErrorCode ierr;
1957: PetscNewLog(fvm,&b);
1958: fvm->data = b;
1960: PetscFVInitialize_Upwind(fvm);
1961: return(0);
1962: }
1964: #include <petscblaslapack.h>
1966: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1967: {
1968: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1969: PetscErrorCode ierr;
1972: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1973: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1974: PetscFree(ls);
1975: return(0);
1976: }
1978: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1979: {
1980: PetscInt Nc, c;
1981: PetscViewerFormat format;
1982: PetscErrorCode ierr;
1985: PetscFVGetNumComponents(fv, &Nc);
1986: PetscViewerGetFormat(viewer, &format);
1987: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1988: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1989: for (c = 0; c < Nc; c++) {
1990: if (fv->componentNames[c]) {
1991: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1992: }
1993: }
1994: return(0);
1995: }
1997: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1998: {
1999: PetscBool iascii;
2005: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2006: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2007: return(0);
2008: }
2010: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2011: {
2013: return(0);
2014: }
2016: /* Overwrites A. Can only handle full-rank problems with m>=n */
2017: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2018: {
2019: PetscBool debug = PETSC_FALSE;
2021: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2022: PetscScalar *R,*Q,*Aback,Alpha;
2025: if (debug) {
2026: PetscMalloc1(m*n,&Aback);
2027: PetscArraycpy(Aback,A,m*n);
2028: }
2030: PetscBLASIntCast(m,&M);
2031: PetscBLASIntCast(n,&N);
2032: PetscBLASIntCast(mstride,&lda);
2033: PetscBLASIntCast(worksize,&ldwork);
2034: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2035: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2036: PetscFPTrapPop();
2037: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2038: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2040: /* Extract an explicit representation of Q */
2041: Q = Ainv;
2042: PetscArraycpy(Q,A,mstride*n);
2043: K = N; /* full rank */
2044: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2045: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2047: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2048: Alpha = 1.0;
2049: ldb = lda;
2050: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2051: /* Ainv is Q, overwritten with inverse */
2053: if (debug) { /* Check that pseudo-inverse worked */
2054: PetscScalar Beta = 0.0;
2055: PetscBLASInt ldc;
2056: K = N;
2057: ldc = N;
2058: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2059: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2060: PetscFree(Aback);
2061: }
2062: return(0);
2063: }
2065: /* Overwrites A. Can handle degenerate problems and m<n. */
2066: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2067: {
2068: PetscBool debug = PETSC_FALSE;
2069: PetscScalar *Brhs, *Aback;
2070: PetscScalar *tmpwork;
2071: PetscReal rcond;
2072: #if defined (PETSC_USE_COMPLEX)
2073: PetscInt rworkSize;
2074: PetscReal *rwork;
2075: #endif
2076: PetscInt i, j, maxmn;
2077: PetscBLASInt M, N, lda, ldb, ldwork;
2078: PetscBLASInt nrhs, irank, info;
2082: if (debug) {
2083: PetscMalloc1(m*n,&Aback);
2084: PetscArraycpy(Aback,A,m*n);
2085: }
2087: /* initialize to identity */
2088: tmpwork = work;
2089: Brhs = Ainv;
2090: maxmn = PetscMax(m,n);
2091: for (j=0; j<maxmn; j++) {
2092: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2093: }
2095: PetscBLASIntCast(m,&M);
2096: PetscBLASIntCast(n,&N);
2097: PetscBLASIntCast(mstride,&lda);
2098: PetscBLASIntCast(maxmn,&ldb);
2099: PetscBLASIntCast(worksize,&ldwork);
2100: rcond = -1;
2101: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2102: nrhs = M;
2103: #if defined(PETSC_USE_COMPLEX)
2104: rworkSize = 5 * PetscMin(M,N);
2105: PetscMalloc1(rworkSize,&rwork);
2106: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info));
2107: PetscFPTrapPop();
2108: PetscFree(rwork);
2109: #else
2110: nrhs = M;
2111: PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info));
2112: PetscFPTrapPop();
2113: #endif
2114: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2115: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2116: 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");
2117: return(0);
2118: }
2120: #if 0
2121: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2122: {
2123: PetscReal grad[2] = {0, 0};
2124: const PetscInt *faces;
2125: PetscInt numFaces, f;
2126: PetscErrorCode ierr;
2129: DMPlexGetConeSize(dm, cell, &numFaces);
2130: DMPlexGetCone(dm, cell, &faces);
2131: for (f = 0; f < numFaces; ++f) {
2132: const PetscInt *fcells;
2133: const CellGeom *cg1;
2134: const FaceGeom *fg;
2136: DMPlexGetSupport(dm, faces[f], &fcells);
2137: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2138: for (i = 0; i < 2; ++i) {
2139: PetscScalar du;
2141: if (fcells[i] == c) continue;
2142: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2143: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2144: grad[0] += fg->grad[!i][0] * du;
2145: grad[1] += fg->grad[!i][1] * du;
2146: }
2147: }
2148: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2149: return(0);
2150: }
2151: #endif
2153: /*
2154: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2156: Input Parameters:
2157: + fvm - The PetscFV object
2158: . numFaces - The number of cell faces which are not constrained
2159: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2161: Level: developer
2163: .seealso: PetscFVCreate()
2164: */
2165: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2166: {
2167: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2168: const PetscBool useSVD = PETSC_TRUE;
2169: const PetscInt maxFaces = ls->maxFaces;
2170: PetscInt dim, f, d;
2171: PetscErrorCode ierr;
2174: if (numFaces > maxFaces) {
2175: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2176: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2177: }
2178: PetscFVGetSpatialDimension(fvm, &dim);
2179: for (f = 0; f < numFaces; ++f) {
2180: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2181: }
2182: /* Overwrites B with garbage, returns Binv in row-major format */
2183: if (useSVD) {
2184: PetscInt maxmn = PetscMax(numFaces, dim);
2185: PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2186: /* Binv shaped in column-major, coldim=maxmn.*/
2187: for (f = 0; f < numFaces; ++f) {
2188: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d + maxmn*f];
2189: }
2190: } else {
2191: PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2192: /* Binv shaped in row-major, rowdim=maxFaces.*/
2193: for (f = 0; f < numFaces; ++f) {
2194: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces + f];
2195: }
2196: }
2197: return(0);
2198: }
2200: /*
2201: neighborVol[f*2+0] contains the left geom
2202: neighborVol[f*2+1] contains the right geom
2203: */
2204: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2205: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2206: {
2207: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2208: void *rctx;
2209: PetscScalar *flux = fvm->fluxWork;
2210: const PetscScalar *constants;
2211: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2212: PetscErrorCode ierr;
2215: PetscDSGetTotalComponents(prob, &Nc);
2216: PetscDSGetTotalDimension(prob, &totDim);
2217: PetscDSGetFieldOffset(prob, field, &off);
2218: PetscDSGetRiemannSolver(prob, field, &riemann);
2219: PetscDSGetContext(prob, field, &rctx);
2220: PetscDSGetConstants(prob, &numConstants, &constants);
2221: PetscFVGetSpatialDimension(fvm, &dim);
2222: PetscFVGetNumComponents(fvm, &pdim);
2223: for (f = 0; f < Nf; ++f) {
2224: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2225: for (d = 0; d < pdim; ++d) {
2226: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2227: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2228: }
2229: }
2230: return(0);
2231: }
2233: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2234: {
2235: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2236: PetscInt dim,m,n,nrhs,minmn,maxmn;
2237: PetscErrorCode ierr;
2241: PetscFVGetSpatialDimension(fvm, &dim);
2242: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2243: ls->maxFaces = maxFaces;
2244: m = ls->maxFaces;
2245: n = dim;
2246: nrhs = ls->maxFaces;
2247: minmn = PetscMin(m,n);
2248: maxmn = PetscMax(m,n);
2249: ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2250: PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);
2251: return(0);
2252: }
2254: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2255: {
2257: fvm->ops->setfromoptions = NULL;
2258: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2259: fvm->ops->view = PetscFVView_LeastSquares;
2260: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2261: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2262: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2263: return(0);
2264: }
2266: /*MC
2267: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2269: Level: intermediate
2271: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2272: M*/
2274: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2275: {
2276: PetscFV_LeastSquares *ls;
2277: PetscErrorCode ierr;
2281: PetscNewLog(fvm, &ls);
2282: fvm->data = ls;
2284: ls->maxFaces = -1;
2285: ls->workSize = -1;
2286: ls->B = NULL;
2287: ls->Binv = NULL;
2288: ls->tau = NULL;
2289: ls->work = NULL;
2291: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2292: PetscFVInitialize_LeastSquares(fvm);
2293: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2294: return(0);
2295: }
2297: /*@
2298: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2300: Not collective
2302: Input parameters:
2303: + fvm - The PetscFV object
2304: - maxFaces - The maximum number of cell faces
2306: Level: intermediate
2308: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2309: @*/
2310: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2311: {
2316: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2317: return(0);
2318: }