Actual source code: fv.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/petscfvimpl.h>
2: #include <petsc/private/dmpleximpl.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: PetscLimiterView - Views a PetscLimiter
129: Collective on lim
131: Input Parameter:
132: + lim - the PetscLimiter object to view
133: - v - the viewer
135: Level: beginner
137: .seealso: PetscLimiterDestroy()
138: @*/
139: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
140: {
145: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
146: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
147: return(0);
148: }
150: /*@
151: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
153: Collective on lim
155: Input Parameter:
156: . lim - the PetscLimiter object to set options for
158: Level: intermediate
160: .seealso: PetscLimiterView()
161: @*/
162: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
163: {
164: const char *defaultType;
165: char name[256];
166: PetscBool flg;
171: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
172: else defaultType = ((PetscObject) lim)->type_name;
173: PetscLimiterRegisterAll();
175: PetscObjectOptionsBegin((PetscObject) lim);
176: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
177: if (flg) {
178: PetscLimiterSetType(lim, name);
179: } else if (!((PetscObject) lim)->type_name) {
180: PetscLimiterSetType(lim, defaultType);
181: }
182: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
183: /* process any options handlers added with PetscObjectAddOptionsHandler() */
184: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
185: PetscOptionsEnd();
186: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
187: return(0);
188: }
190: /*@C
191: PetscLimiterSetUp - Construct data structures for the PetscLimiter
193: Collective on lim
195: Input Parameter:
196: . lim - the PetscLimiter object to setup
198: Level: intermediate
200: .seealso: PetscLimiterView(), PetscLimiterDestroy()
201: @*/
202: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
203: {
208: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
209: return(0);
210: }
212: /*@
213: PetscLimiterDestroy - Destroys a PetscLimiter object
215: Collective on lim
217: Input Parameter:
218: . lim - the PetscLimiter object to destroy
220: Level: beginner
222: .seealso: PetscLimiterView()
223: @*/
224: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
225: {
229: if (!*lim) return(0);
232: if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
233: ((PetscObject) (*lim))->refct = 0;
235: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
236: PetscHeaderDestroy(lim);
237: return(0);
238: }
240: /*@
241: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
243: Collective
245: Input Parameter:
246: . comm - The communicator for the PetscLimiter object
248: Output Parameter:
249: . lim - The PetscLimiter object
251: Level: beginner
253: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
254: @*/
255: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
256: {
257: PetscLimiter l;
262: PetscCitationsRegister(LimiterCitation,&Limitercite);
263: *lim = NULL;
264: PetscFVInitializePackage();
266: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
268: *lim = l;
269: return(0);
270: }
272: /*@
273: PetscLimiterLimit - Limit the flux
275: Input Parameters:
276: + lim - The PetscLimiter
277: - flim - The input field
279: Output Parameter:
280: . phi - The limited field
282: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
283: $ The classical flux-limited formulation is psi(r) where
284: $
285: $ r = (u[0] - u[-1]) / (u[1] - u[0])
286: $
287: $ The second order TVD region is bounded by
288: $
289: $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
290: $
291: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
292: $ phi(r)(r+1)/2 in which the reconstructed interface values are
293: $
294: $ u(v) = u[0] + phi(r) (grad u)[0] v
295: $
296: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
297: $
298: $ 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))
299: $
300: $ For a nicer symmetric formulation, rewrite in terms of
301: $
302: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
303: $
304: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
305: $
306: $ phi(r) = phi(1/r)
307: $
308: $ becomes
309: $
310: $ w(f) = w(1-f).
311: $
312: $ The limiters below implement this final form w(f). The reference methods are
313: $
314: $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
316: Level: beginner
318: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
319: @*/
320: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
321: {
327: (*lim->ops->limit)(lim, flim, phi);
328: return(0);
329: }
331: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
332: {
333: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
334: PetscErrorCode ierr;
337: PetscFree(l);
338: return(0);
339: }
341: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
342: {
343: PetscViewerFormat format;
344: PetscErrorCode ierr;
347: PetscViewerGetFormat(viewer, &format);
348: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
349: return(0);
350: }
352: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
353: {
354: PetscBool iascii;
360: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
361: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
362: return(0);
363: }
365: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
366: {
368: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
369: return(0);
370: }
372: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
373: {
375: lim->ops->view = PetscLimiterView_Sin;
376: lim->ops->destroy = PetscLimiterDestroy_Sin;
377: lim->ops->limit = PetscLimiterLimit_Sin;
378: return(0);
379: }
381: /*MC
382: PETSCLIMITERSIN = "sin" - A PetscLimiter object
384: Level: intermediate
386: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
387: M*/
389: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
390: {
391: PetscLimiter_Sin *l;
392: PetscErrorCode ierr;
396: PetscNewLog(lim, &l);
397: lim->data = l;
399: PetscLimiterInitialize_Sin(lim);
400: return(0);
401: }
403: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
404: {
405: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
406: PetscErrorCode ierr;
409: PetscFree(l);
410: return(0);
411: }
413: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
414: {
415: PetscViewerFormat format;
416: PetscErrorCode ierr;
419: PetscViewerGetFormat(viewer, &format);
420: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
421: return(0);
422: }
424: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
425: {
426: PetscBool iascii;
432: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
433: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
434: return(0);
435: }
437: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
438: {
440: *phi = 0.0;
441: return(0);
442: }
444: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
445: {
447: lim->ops->view = PetscLimiterView_Zero;
448: lim->ops->destroy = PetscLimiterDestroy_Zero;
449: lim->ops->limit = PetscLimiterLimit_Zero;
450: return(0);
451: }
453: /*MC
454: PETSCLIMITERZERO = "zero" - A PetscLimiter object
456: Level: intermediate
458: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
459: M*/
461: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
462: {
463: PetscLimiter_Zero *l;
464: PetscErrorCode ierr;
468: PetscNewLog(lim, &l);
469: lim->data = l;
471: PetscLimiterInitialize_Zero(lim);
472: return(0);
473: }
475: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
476: {
477: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
478: PetscErrorCode ierr;
481: PetscFree(l);
482: return(0);
483: }
485: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
486: {
487: PetscViewerFormat format;
488: PetscErrorCode ierr;
491: PetscViewerGetFormat(viewer, &format);
492: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
493: return(0);
494: }
496: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
497: {
498: PetscBool iascii;
504: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
505: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
506: return(0);
507: }
509: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
510: {
512: *phi = 1.0;
513: return(0);
514: }
516: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
517: {
519: lim->ops->view = PetscLimiterView_None;
520: lim->ops->destroy = PetscLimiterDestroy_None;
521: lim->ops->limit = PetscLimiterLimit_None;
522: return(0);
523: }
525: /*MC
526: PETSCLIMITERNONE = "none" - A PetscLimiter object
528: Level: intermediate
530: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
531: M*/
533: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
534: {
535: PetscLimiter_None *l;
536: PetscErrorCode ierr;
540: PetscNewLog(lim, &l);
541: lim->data = l;
543: PetscLimiterInitialize_None(lim);
544: return(0);
545: }
547: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
548: {
549: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
550: PetscErrorCode ierr;
553: PetscFree(l);
554: return(0);
555: }
557: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
558: {
559: PetscViewerFormat format;
560: PetscErrorCode ierr;
563: PetscViewerGetFormat(viewer, &format);
564: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
565: return(0);
566: }
568: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
569: {
570: PetscBool iascii;
576: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
577: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
578: return(0);
579: }
581: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
582: {
584: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
585: return(0);
586: }
588: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
589: {
591: lim->ops->view = PetscLimiterView_Minmod;
592: lim->ops->destroy = PetscLimiterDestroy_Minmod;
593: lim->ops->limit = PetscLimiterLimit_Minmod;
594: return(0);
595: }
597: /*MC
598: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
600: Level: intermediate
602: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
603: M*/
605: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
606: {
607: PetscLimiter_Minmod *l;
608: PetscErrorCode ierr;
612: PetscNewLog(lim, &l);
613: lim->data = l;
615: PetscLimiterInitialize_Minmod(lim);
616: return(0);
617: }
619: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
620: {
621: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
622: PetscErrorCode ierr;
625: PetscFree(l);
626: return(0);
627: }
629: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
630: {
631: PetscViewerFormat format;
632: PetscErrorCode ierr;
635: PetscViewerGetFormat(viewer, &format);
636: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
637: return(0);
638: }
640: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
641: {
642: PetscBool iascii;
648: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
649: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
650: return(0);
651: }
653: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
654: {
656: *phi = PetscMax(0, 4*f*(1-f));
657: return(0);
658: }
660: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
661: {
663: lim->ops->view = PetscLimiterView_VanLeer;
664: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
665: lim->ops->limit = PetscLimiterLimit_VanLeer;
666: return(0);
667: }
669: /*MC
670: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
672: Level: intermediate
674: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
675: M*/
677: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
678: {
679: PetscLimiter_VanLeer *l;
680: PetscErrorCode ierr;
684: PetscNewLog(lim, &l);
685: lim->data = l;
687: PetscLimiterInitialize_VanLeer(lim);
688: return(0);
689: }
691: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
692: {
693: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
694: PetscErrorCode ierr;
697: PetscFree(l);
698: return(0);
699: }
701: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
702: {
703: PetscViewerFormat format;
704: PetscErrorCode ierr;
707: PetscViewerGetFormat(viewer, &format);
708: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
709: return(0);
710: }
712: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
713: {
714: PetscBool iascii;
720: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
721: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
722: return(0);
723: }
725: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
726: {
728: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
729: return(0);
730: }
732: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
733: {
735: lim->ops->view = PetscLimiterView_VanAlbada;
736: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
737: lim->ops->limit = PetscLimiterLimit_VanAlbada;
738: return(0);
739: }
741: /*MC
742: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
744: Level: intermediate
746: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
747: M*/
749: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
750: {
751: PetscLimiter_VanAlbada *l;
752: PetscErrorCode ierr;
756: PetscNewLog(lim, &l);
757: lim->data = l;
759: PetscLimiterInitialize_VanAlbada(lim);
760: return(0);
761: }
763: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
764: {
765: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
766: PetscErrorCode ierr;
769: PetscFree(l);
770: return(0);
771: }
773: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
774: {
775: PetscViewerFormat format;
776: PetscErrorCode ierr;
779: PetscViewerGetFormat(viewer, &format);
780: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
781: return(0);
782: }
784: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
785: {
786: PetscBool iascii;
792: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
793: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
794: return(0);
795: }
797: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
798: {
800: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
801: return(0);
802: }
804: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
805: {
807: lim->ops->view = PetscLimiterView_Superbee;
808: lim->ops->destroy = PetscLimiterDestroy_Superbee;
809: lim->ops->limit = PetscLimiterLimit_Superbee;
810: return(0);
811: }
813: /*MC
814: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
816: Level: intermediate
818: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
819: M*/
821: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
822: {
823: PetscLimiter_Superbee *l;
824: PetscErrorCode ierr;
828: PetscNewLog(lim, &l);
829: lim->data = l;
831: PetscLimiterInitialize_Superbee(lim);
832: return(0);
833: }
835: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
836: {
837: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
838: PetscErrorCode ierr;
841: PetscFree(l);
842: return(0);
843: }
845: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
846: {
847: PetscViewerFormat format;
848: PetscErrorCode ierr;
851: PetscViewerGetFormat(viewer, &format);
852: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
853: return(0);
854: }
856: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
857: {
858: PetscBool iascii;
864: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
865: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
866: return(0);
867: }
869: /* aka Barth-Jespersen */
870: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
871: {
873: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
874: return(0);
875: }
877: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
878: {
880: lim->ops->view = PetscLimiterView_MC;
881: lim->ops->destroy = PetscLimiterDestroy_MC;
882: lim->ops->limit = PetscLimiterLimit_MC;
883: return(0);
884: }
886: /*MC
887: PETSCLIMITERMC = "mc" - A PetscLimiter object
889: Level: intermediate
891: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
892: M*/
894: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
895: {
896: PetscLimiter_MC *l;
897: PetscErrorCode ierr;
901: PetscNewLog(lim, &l);
902: lim->data = l;
904: PetscLimiterInitialize_MC(lim);
905: return(0);
906: }
908: PetscClassId PETSCFV_CLASSID = 0;
910: PetscFunctionList PetscFVList = NULL;
911: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
913: /*@C
914: PetscFVRegister - Adds a new PetscFV implementation
916: Not Collective
918: Input Parameters:
919: + name - The name of a new user-defined creation routine
920: - create_func - The creation routine itself
922: Notes:
923: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
925: Sample usage:
926: .vb
927: PetscFVRegister("my_fv", MyPetscFVCreate);
928: .ve
930: Then, your PetscFV type can be chosen with the procedural interface via
931: .vb
932: PetscFVCreate(MPI_Comm, PetscFV *);
933: PetscFVSetType(PetscFV, "my_fv");
934: .ve
935: or at runtime via the option
936: .vb
937: -petscfv_type my_fv
938: .ve
940: Level: advanced
942: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
944: @*/
945: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
946: {
950: PetscFunctionListAdd(&PetscFVList, sname, function);
951: return(0);
952: }
954: /*@C
955: PetscFVSetType - Builds a particular PetscFV
957: Collective on fvm
959: Input Parameters:
960: + fvm - The PetscFV object
961: - name - The kind of FVM space
963: Options Database Key:
964: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
966: Level: intermediate
968: .seealso: PetscFVGetType(), PetscFVCreate()
969: @*/
970: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
971: {
972: PetscErrorCode (*r)(PetscFV);
973: PetscBool match;
978: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
979: if (match) return(0);
981: PetscFVRegisterAll();
982: PetscFunctionListFind(PetscFVList, name, &r);
983: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
985: if (fvm->ops->destroy) {
986: (*fvm->ops->destroy)(fvm);
987: fvm->ops->destroy = NULL;
988: }
989: (*r)(fvm);
990: PetscObjectChangeTypeName((PetscObject) fvm, name);
991: return(0);
992: }
994: /*@C
995: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
997: Not Collective
999: Input Parameter:
1000: . fvm - The PetscFV
1002: Output Parameter:
1003: . name - The PetscFV type name
1005: Level: intermediate
1007: .seealso: PetscFVSetType(), PetscFVCreate()
1008: @*/
1009: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1010: {
1016: PetscFVRegisterAll();
1017: *name = ((PetscObject) fvm)->type_name;
1018: return(0);
1019: }
1021: /*@C
1022: PetscFVView - Views a PetscFV
1024: Collective on fvm
1026: Input Parameter:
1027: + fvm - the PetscFV object to view
1028: - v - the viewer
1030: Level: beginner
1032: .seealso: PetscFVDestroy()
1033: @*/
1034: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1035: {
1040: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1041: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1042: return(0);
1043: }
1045: /*@
1046: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1048: Collective on fvm
1050: Input Parameter:
1051: . fvm - the PetscFV object to set options for
1053: Options Database Key:
1054: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1056: Level: intermediate
1058: .seealso: PetscFVView()
1059: @*/
1060: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1061: {
1062: const char *defaultType;
1063: char name[256];
1064: PetscBool flg;
1069: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1070: else defaultType = ((PetscObject) fvm)->type_name;
1071: PetscFVRegisterAll();
1073: PetscObjectOptionsBegin((PetscObject) fvm);
1074: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1075: if (flg) {
1076: PetscFVSetType(fvm, name);
1077: } else if (!((PetscObject) fvm)->type_name) {
1078: PetscFVSetType(fvm, defaultType);
1080: }
1081: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1082: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1083: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1084: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1085: PetscLimiterSetFromOptions(fvm->limiter);
1086: PetscOptionsEnd();
1087: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1088: return(0);
1089: }
1091: /*@
1092: PetscFVSetUp - Construct data structures for the PetscFV
1094: Collective on fvm
1096: Input Parameter:
1097: . fvm - the PetscFV object to setup
1099: Level: intermediate
1101: .seealso: PetscFVView(), PetscFVDestroy()
1102: @*/
1103: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1104: {
1109: PetscLimiterSetUp(fvm->limiter);
1110: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1111: return(0);
1112: }
1114: /*@
1115: PetscFVDestroy - Destroys a PetscFV object
1117: Collective on fvm
1119: Input Parameter:
1120: . fvm - the PetscFV object to destroy
1122: Level: beginner
1124: .seealso: PetscFVView()
1125: @*/
1126: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1127: {
1128: PetscInt i;
1132: if (!*fvm) return(0);
1135: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1136: ((PetscObject) (*fvm))->refct = 0;
1138: for (i = 0; i < (*fvm)->numComponents; i++) {
1139: PetscFree((*fvm)->componentNames[i]);
1140: }
1141: PetscFree((*fvm)->componentNames);
1142: PetscLimiterDestroy(&(*fvm)->limiter);
1143: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1144: PetscFree((*fvm)->fluxWork);
1145: PetscQuadratureDestroy(&(*fvm)->quadrature);
1146: PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);
1148: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1149: PetscHeaderDestroy(fvm);
1150: return(0);
1151: }
1153: /*@
1154: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1156: Collective
1158: Input Parameter:
1159: . comm - The communicator for the PetscFV object
1161: Output Parameter:
1162: . fvm - The PetscFV object
1164: Level: beginner
1166: .seealso: PetscFVSetType(), PETSCFVUPWIND
1167: @*/
1168: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1169: {
1170: PetscFV f;
1175: *fvm = NULL;
1176: PetscFVInitializePackage();
1178: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1179: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1181: PetscLimiterCreate(comm, &f->limiter);
1182: f->numComponents = 1;
1183: f->dim = 0;
1184: f->computeGradients = PETSC_FALSE;
1185: f->fluxWork = NULL;
1186: PetscCalloc1(f->numComponents,&f->componentNames);
1188: *fvm = f;
1189: return(0);
1190: }
1192: /*@
1193: PetscFVSetLimiter - Set the limiter object
1195: Logically collective on fvm
1197: Input Parameters:
1198: + fvm - the PetscFV object
1199: - lim - The PetscLimiter
1201: Level: intermediate
1203: .seealso: PetscFVGetLimiter()
1204: @*/
1205: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1206: {
1212: PetscLimiterDestroy(&fvm->limiter);
1213: PetscObjectReference((PetscObject) lim);
1214: fvm->limiter = lim;
1215: return(0);
1216: }
1218: /*@
1219: PetscFVGetLimiter - Get the limiter object
1221: Not collective
1223: Input Parameter:
1224: . fvm - the PetscFV object
1226: Output Parameter:
1227: . lim - The PetscLimiter
1229: Level: intermediate
1231: .seealso: PetscFVSetLimiter()
1232: @*/
1233: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1234: {
1238: *lim = fvm->limiter;
1239: return(0);
1240: }
1242: /*@
1243: PetscFVSetNumComponents - Set the number of field components
1245: Logically collective on fvm
1247: Input Parameters:
1248: + fvm - the PetscFV object
1249: - comp - The number of components
1251: Level: intermediate
1253: .seealso: PetscFVGetNumComponents()
1254: @*/
1255: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1256: {
1261: if (fvm->numComponents != comp) {
1262: PetscInt i;
1264: for (i = 0; i < fvm->numComponents; i++) {
1265: PetscFree(fvm->componentNames[i]);
1266: }
1267: PetscFree(fvm->componentNames);
1268: PetscCalloc1(comp,&fvm->componentNames);
1269: }
1270: fvm->numComponents = comp;
1271: PetscFree(fvm->fluxWork);
1272: PetscMalloc1(comp, &fvm->fluxWork);
1273: return(0);
1274: }
1276: /*@
1277: PetscFVGetNumComponents - Get the number of field components
1279: Not collective
1281: Input Parameter:
1282: . fvm - the PetscFV object
1284: Output Parameter:
1285: , comp - The number of components
1287: Level: intermediate
1289: .seealso: PetscFVSetNumComponents()
1290: @*/
1291: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1292: {
1296: *comp = fvm->numComponents;
1297: return(0);
1298: }
1300: /*@C
1301: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1303: Logically collective on fvm
1304: Input Parameters:
1305: + fvm - the PetscFV object
1306: . comp - the component number
1307: - name - the component name
1309: Level: intermediate
1311: .seealso: PetscFVGetComponentName()
1312: @*/
1313: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1314: {
1318: PetscFree(fvm->componentNames[comp]);
1319: PetscStrallocpy(name,&fvm->componentNames[comp]);
1320: return(0);
1321: }
1323: /*@C
1324: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1326: Logically collective on fvm
1327: Input Parameters:
1328: + fvm - the PetscFV object
1329: - comp - the component number
1331: Output Parameter:
1332: . name - the component name
1334: Level: intermediate
1336: .seealso: PetscFVSetComponentName()
1337: @*/
1338: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1339: {
1341: *name = fvm->componentNames[comp];
1342: return(0);
1343: }
1345: /*@
1346: PetscFVSetSpatialDimension - Set the spatial dimension
1348: Logically collective on fvm
1350: Input Parameters:
1351: + fvm - the PetscFV object
1352: - dim - The spatial dimension
1354: Level: intermediate
1356: .seealso: PetscFVGetSpatialDimension()
1357: @*/
1358: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1359: {
1362: fvm->dim = dim;
1363: return(0);
1364: }
1366: /*@
1367: PetscFVGetSpatialDimension - Get the spatial dimension
1369: Logically collective on fvm
1371: Input Parameter:
1372: . fvm - the PetscFV object
1374: Output Parameter:
1375: . dim - The spatial dimension
1377: Level: intermediate
1379: .seealso: PetscFVSetSpatialDimension()
1380: @*/
1381: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1382: {
1386: *dim = fvm->dim;
1387: return(0);
1388: }
1390: /*@
1391: PetscFVSetComputeGradients - Toggle computation of cell gradients
1393: Logically collective on fvm
1395: Input Parameters:
1396: + fvm - the PetscFV object
1397: - computeGradients - Flag to compute cell gradients
1399: Level: intermediate
1401: .seealso: PetscFVGetComputeGradients()
1402: @*/
1403: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1404: {
1407: fvm->computeGradients = computeGradients;
1408: return(0);
1409: }
1411: /*@
1412: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1414: Not collective
1416: Input Parameter:
1417: . fvm - the PetscFV object
1419: Output Parameter:
1420: . computeGradients - Flag to compute cell gradients
1422: Level: intermediate
1424: .seealso: PetscFVSetComputeGradients()
1425: @*/
1426: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1427: {
1431: *computeGradients = fvm->computeGradients;
1432: return(0);
1433: }
1435: /*@
1436: PetscFVSetQuadrature - Set the quadrature object
1438: Logically collective on fvm
1440: Input Parameters:
1441: + fvm - the PetscFV object
1442: - q - The PetscQuadrature
1444: Level: intermediate
1446: .seealso: PetscFVGetQuadrature()
1447: @*/
1448: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1449: {
1454: PetscQuadratureDestroy(&fvm->quadrature);
1455: PetscObjectReference((PetscObject) q);
1456: fvm->quadrature = q;
1457: return(0);
1458: }
1460: /*@
1461: PetscFVGetQuadrature - Get the quadrature object
1463: Not collective
1465: Input Parameter:
1466: . fvm - the PetscFV object
1468: Output Parameter:
1469: . lim - The PetscQuadrature
1471: Level: intermediate
1473: .seealso: PetscFVSetQuadrature()
1474: @*/
1475: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1476: {
1480: if (!fvm->quadrature) {
1481: /* Create default 1-point quadrature */
1482: PetscReal *points, *weights;
1485: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1486: PetscCalloc1(fvm->dim, &points);
1487: PetscMalloc1(1, &weights);
1488: weights[0] = 1.0;
1489: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1490: }
1491: *q = fvm->quadrature;
1492: return(0);
1493: }
1495: /*@
1496: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1498: Not collective
1500: Input Parameter:
1501: . fvm - The PetscFV object
1503: Output Parameter:
1504: . sp - The PetscDualSpace object
1506: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1508: Level: intermediate
1510: .seealso: PetscFVCreate()
1511: @*/
1512: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1513: {
1517: if (!fvm->dualSpace) {
1518: DM K;
1519: PetscInt dim, Nc, c;
1520: PetscErrorCode ierr;
1522: PetscFVGetSpatialDimension(fvm, &dim);
1523: PetscFVGetNumComponents(fvm, &Nc);
1524: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1525: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1526: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1527: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1528: PetscDualSpaceSetDM(fvm->dualSpace, K);
1529: DMDestroy(&K);
1530: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1531: /* Should we be using PetscFVGetQuadrature() here? */
1532: for (c = 0; c < Nc; ++c) {
1533: PetscQuadrature qc;
1534: PetscReal *points, *weights;
1535: PetscErrorCode ierr;
1537: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1538: PetscCalloc1(dim, &points);
1539: PetscCalloc1(Nc, &weights);
1540: weights[c] = 1.0;
1541: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1542: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1543: PetscQuadratureDestroy(&qc);
1544: }
1545: }
1546: *sp = fvm->dualSpace;
1547: return(0);
1548: }
1550: /*@
1551: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1553: Not collective
1555: Input Parameters:
1556: + fvm - The PetscFV object
1557: - sp - The PetscDualSpace object
1559: Level: intermediate
1561: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1563: .seealso: PetscFVCreate()
1564: @*/
1565: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1566: {
1572: PetscDualSpaceDestroy(&fvm->dualSpace);
1573: fvm->dualSpace = sp;
1574: PetscObjectReference((PetscObject) fvm->dualSpace);
1575: return(0);
1576: }
1578: /*@C
1579: PetscFVGetDefaultTabulation - Returns the tabulation of the basis functions at the quadrature points
1581: Not collective
1583: Input Parameter:
1584: . fvm - The PetscFV object
1586: Output Parameters:
1587: + B - The basis function values at quadrature points
1588: . D - The basis function derivatives at quadrature points
1589: - H - The basis function second derivatives at quadrature points
1591: Note:
1592: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1593: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1594: $ 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
1596: Level: intermediate
1598: .seealso: PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1599: @*/
1600: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1601: {
1602: PetscInt npoints;
1603: const PetscReal *points;
1604: PetscErrorCode ierr;
1611: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1612: if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1613: if (B) *B = fvm->B;
1614: if (D) *D = fvm->D;
1615: if (H) *H = fvm->H;
1616: return(0);
1617: }
1619: /*@C
1620: PetscFVGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1622: Not collective
1624: Input Parameters:
1625: + fvm - The PetscFV object
1626: . npoints - The number of tabulation points
1627: - points - The tabulation point coordinates
1629: Output Parameters:
1630: + B - The basis function values at tabulation points
1631: . D - The basis function derivatives at tabulation points
1632: - H - The basis function second derivatives at tabulation points
1634: Note:
1635: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1636: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1637: $ 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
1639: Level: intermediate
1641: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
1642: @*/
1643: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1644: {
1645: PetscInt pdim = 1; /* Dimension of approximation space P */
1646: PetscInt dim; /* Spatial dimension */
1647: PetscInt comp; /* Field components */
1648: PetscInt p, d, c, e;
1649: PetscErrorCode ierr;
1657: PetscFVGetSpatialDimension(fvm, &dim);
1658: PetscFVGetNumComponents(fvm, &comp);
1659: if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1660: if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1661: if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1662: if (B) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) (*B)[(p*pdim + d)*comp + c] = 1.0;}
1663: if (D) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim; ++e) (*D)[((p*pdim + d)*comp + c)*dim + e] = 0.0;}
1664: if (H) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim*dim; ++e) (*H)[((p*pdim + d)*comp + c)*dim*dim + e] = 0.0;}
1665: return(0);
1666: }
1668: /*@C
1669: PetscFVRestoreTabulation - Frees memory from the associated tabulation.
1671: Not collective
1673: Input Parameters:
1674: + fvm - The PetscFV object
1675: . npoints - The number of tabulation points
1676: . points - The tabulation point coordinates
1677: . B - The basis function values at tabulation points
1678: . D - The basis function derivatives at tabulation points
1679: - H - The basis function second derivatives at tabulation points
1681: Note:
1682: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1683: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1684: $ 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
1686: Level: intermediate
1688: .seealso: PetscFVGetTabulation(), PetscFVGetDefaultTabulation()
1689: @*/
1690: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1691: {
1696: if (B && *B) {PetscFree(*B);}
1697: if (D && *D) {PetscFree(*D);}
1698: if (H && *H) {PetscFree(*H);}
1699: return(0);
1700: }
1702: /*@C
1703: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1705: Input Parameters:
1706: + fvm - The PetscFV object
1707: . numFaces - The number of cell faces which are not constrained
1708: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1710: Level: advanced
1712: .seealso: PetscFVCreate()
1713: @*/
1714: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1715: {
1720: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1721: return(0);
1722: }
1724: /*@C
1725: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1727: Not collective
1729: Input Parameters:
1730: + fvm - The PetscFV object for the field being integrated
1731: . prob - The PetscDS specifing the discretizations and continuum functions
1732: . field - The field being integrated
1733: . Nf - The number of faces in the chunk
1734: . fgeom - The face geometry for each face in the chunk
1735: . neighborVol - The volume for each pair of cells in the chunk
1736: . uL - The state from the cell on the left
1737: - uR - The state from the cell on the right
1739: Output Parameter
1740: + fluxL - the left fluxes for each face
1741: - fluxR - the right fluxes for each face
1743: Level: developer
1745: .seealso: PetscFVCreate()
1746: @*/
1747: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1748: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1749: {
1754: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1755: return(0);
1756: }
1758: /*@
1759: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1760: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1761: sparsity). It is also used to create an interpolation between regularly refined meshes.
1763: Input Parameter:
1764: . fv - The initial PetscFV
1766: Output Parameter:
1767: . fvRef - The refined PetscFV
1769: Level: advanced
1771: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1772: @*/
1773: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1774: {
1775: PetscDualSpace Q, Qref;
1776: DM K, Kref;
1777: PetscQuadrature q, qref;
1778: CellRefiner cellRefiner;
1779: PetscReal *v0;
1780: PetscReal *jac, *invjac;
1781: PetscInt numComp, numSubelements, s;
1782: PetscErrorCode ierr;
1785: PetscFVGetDualSpace(fv, &Q);
1786: PetscFVGetQuadrature(fv, &q);
1787: PetscDualSpaceGetDM(Q, &K);
1788: /* Create dual space */
1789: PetscDualSpaceDuplicate(Q, &Qref);
1790: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1791: PetscDualSpaceSetDM(Qref, Kref);
1792: DMDestroy(&Kref);
1793: PetscDualSpaceSetUp(Qref);
1794: /* Create volume */
1795: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1796: PetscFVSetDualSpace(*fvRef, Qref);
1797: PetscFVGetNumComponents(fv, &numComp);
1798: PetscFVSetNumComponents(*fvRef, numComp);
1799: PetscFVSetUp(*fvRef);
1800: /* Create quadrature */
1801: DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1802: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1803: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1804: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1805: for (s = 0; s < numSubelements; ++s) {
1806: PetscQuadrature qs;
1807: const PetscReal *points, *weights;
1808: PetscReal *p, *w;
1809: PetscInt dim, Nc, npoints, np;
1811: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1812: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1813: np = npoints/numSubelements;
1814: PetscMalloc1(np*dim,&p);
1815: PetscMalloc1(np*Nc,&w);
1816: PetscArraycpy(p, &points[s*np*dim], np*dim);
1817: PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1818: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1819: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1820: PetscQuadratureDestroy(&qs);
1821: }
1822: CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1823: PetscFVSetQuadrature(*fvRef, qref);
1824: PetscQuadratureDestroy(&qref);
1825: PetscDualSpaceDestroy(&Qref);
1826: return(0);
1827: }
1829: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1830: {
1831: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1835: PetscFree(b);
1836: return(0);
1837: }
1839: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1840: {
1841: PetscInt Nc, c;
1842: PetscViewerFormat format;
1843: PetscErrorCode ierr;
1846: PetscFVGetNumComponents(fv, &Nc);
1847: PetscViewerGetFormat(viewer, &format);
1848: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1849: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1850: for (c = 0; c < Nc; c++) {
1851: if (fv->componentNames[c]) {
1852: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1853: }
1854: }
1855: return(0);
1856: }
1858: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1859: {
1860: PetscBool iascii;
1866: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1867: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1868: return(0);
1869: }
1871: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1872: {
1874: return(0);
1875: }
1877: /*
1878: neighborVol[f*2+0] contains the left geom
1879: neighborVol[f*2+1] contains the right geom
1880: */
1881: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1882: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1883: {
1884: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1885: void *rctx;
1886: PetscScalar *flux = fvm->fluxWork;
1887: const PetscScalar *constants;
1888: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1889: PetscErrorCode ierr;
1892: PetscDSGetTotalComponents(prob, &Nc);
1893: PetscDSGetTotalDimension(prob, &totDim);
1894: PetscDSGetFieldOffset(prob, field, &off);
1895: PetscDSGetRiemannSolver(prob, field, &riemann);
1896: PetscDSGetContext(prob, field, &rctx);
1897: PetscDSGetConstants(prob, &numConstants, &constants);
1898: PetscFVGetSpatialDimension(fvm, &dim);
1899: PetscFVGetNumComponents(fvm, &pdim);
1900: for (f = 0; f < Nf; ++f) {
1901: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1902: for (d = 0; d < pdim; ++d) {
1903: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1904: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1905: }
1906: }
1907: return(0);
1908: }
1910: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1911: {
1913: fvm->ops->setfromoptions = NULL;
1914: fvm->ops->setup = PetscFVSetUp_Upwind;
1915: fvm->ops->view = PetscFVView_Upwind;
1916: fvm->ops->destroy = PetscFVDestroy_Upwind;
1917: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1918: return(0);
1919: }
1921: /*MC
1922: PETSCFVUPWIND = "upwind" - A PetscFV object
1924: Level: intermediate
1926: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1927: M*/
1929: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1930: {
1931: PetscFV_Upwind *b;
1932: PetscErrorCode ierr;
1936: PetscNewLog(fvm,&b);
1937: fvm->data = b;
1939: PetscFVInitialize_Upwind(fvm);
1940: return(0);
1941: }
1943: #include <petscblaslapack.h>
1945: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1946: {
1947: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1948: PetscErrorCode ierr;
1951: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1952: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1953: PetscFree(ls);
1954: return(0);
1955: }
1957: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1958: {
1959: PetscInt Nc, c;
1960: PetscViewerFormat format;
1961: PetscErrorCode ierr;
1964: PetscFVGetNumComponents(fv, &Nc);
1965: PetscViewerGetFormat(viewer, &format);
1966: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1967: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1968: for (c = 0; c < Nc; c++) {
1969: if (fv->componentNames[c]) {
1970: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1971: }
1972: }
1973: return(0);
1974: }
1976: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1977: {
1978: PetscBool iascii;
1984: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1985: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1986: return(0);
1987: }
1989: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1990: {
1992: return(0);
1993: }
1995: /* Overwrites A. Can only handle full-rank problems with m>=n */
1996: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1997: {
1998: PetscBool debug = PETSC_FALSE;
2000: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2001: PetscScalar *R,*Q,*Aback,Alpha;
2004: if (debug) {
2005: PetscMalloc1(m*n,&Aback);
2006: PetscArraycpy(Aback,A,m*n);
2007: }
2009: PetscBLASIntCast(m,&M);
2010: PetscBLASIntCast(n,&N);
2011: PetscBLASIntCast(mstride,&lda);
2012: PetscBLASIntCast(worksize,&ldwork);
2013: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2014: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2015: PetscFPTrapPop();
2016: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2017: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2019: /* Extract an explicit representation of Q */
2020: Q = Ainv;
2021: PetscArraycpy(Q,A,mstride*n);
2022: K = N; /* full rank */
2023: #if defined(PETSC_MISSING_LAPACK_ORGQR)
2024: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ORGQR - Lapack routine is unavailable.");
2025: #else
2026: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2027: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2028: #endif
2030: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2031: Alpha = 1.0;
2032: ldb = lda;
2033: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2034: /* Ainv is Q, overwritten with inverse */
2036: if (debug) { /* Check that pseudo-inverse worked */
2037: PetscScalar Beta = 0.0;
2038: PetscBLASInt ldc;
2039: K = N;
2040: ldc = N;
2041: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2042: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2043: PetscFree(Aback);
2044: }
2045: return(0);
2046: }
2048: /* Overwrites A. Can handle degenerate problems and m<n. */
2049: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2050: {
2051: PetscBool debug = PETSC_FALSE;
2052: PetscScalar *Brhs, *Aback;
2053: PetscScalar *tmpwork;
2054: PetscReal rcond;
2055: #if defined (PETSC_USE_COMPLEX)
2056: PetscInt rworkSize;
2057: PetscReal *rwork;
2058: #endif
2059: PetscInt i, j, maxmn;
2060: PetscBLASInt M, N, lda, ldb, ldwork;
2061: PetscBLASInt nrhs, irank, info;
2065: if (debug) {
2066: PetscMalloc1(m*n,&Aback);
2067: PetscArraycpy(Aback,A,m*n);
2068: }
2070: /* initialize to identity */
2071: tmpwork = Ainv;
2072: Brhs = work;
2073: maxmn = PetscMax(m,n);
2074: for (j=0; j<maxmn; j++) {
2075: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2076: }
2078: PetscBLASIntCast(m,&M);
2079: PetscBLASIntCast(n,&N);
2080: PetscBLASIntCast(mstride,&lda);
2081: PetscBLASIntCast(maxmn,&ldb);
2082: PetscBLASIntCast(worksize,&ldwork);
2083: rcond = -1;
2084: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2085: nrhs = M;
2086: #if defined(PETSC_USE_COMPLEX)
2087: rworkSize = 5 * PetscMin(M,N);
2088: PetscMalloc1(rworkSize,&rwork);
2089: #if defined(PETSC_MISSING_LAPACK_GELSS)
2090: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GELSS - Lapack routine is unavailable.");
2091: #else
2092: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2093: #endif
2094: PetscFPTrapPop();
2095: PetscFree(rwork);
2096: #else
2097: nrhs = M;
2098: #if defined(PETSC_MISSING_LAPACK_GELSS)
2099: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GELSS - Lapack routine is unavailable.");
2100: #else
2101: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2102: #endif
2103: PetscFPTrapPop();
2104: #endif
2105: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2106: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2107: 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");
2109: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2110: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2111: for (i=0; i<n; i++) {
2112: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2113: }
2114: return(0);
2115: }
2117: #if 0
2118: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2119: {
2120: PetscReal grad[2] = {0, 0};
2121: const PetscInt *faces;
2122: PetscInt numFaces, f;
2123: PetscErrorCode ierr;
2126: DMPlexGetConeSize(dm, cell, &numFaces);
2127: DMPlexGetCone(dm, cell, &faces);
2128: for (f = 0; f < numFaces; ++f) {
2129: const PetscInt *fcells;
2130: const CellGeom *cg1;
2131: const FaceGeom *fg;
2133: DMPlexGetSupport(dm, faces[f], &fcells);
2134: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2135: for (i = 0; i < 2; ++i) {
2136: PetscScalar du;
2138: if (fcells[i] == c) continue;
2139: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2140: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2141: grad[0] += fg->grad[!i][0] * du;
2142: grad[1] += fg->grad[!i][1] * du;
2143: }
2144: }
2145: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2146: return(0);
2147: }
2148: #endif
2150: /*
2151: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2153: Input Parameters:
2154: + fvm - The PetscFV object
2155: . numFaces - The number of cell faces which are not constrained
2156: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2158: Level: developer
2160: .seealso: PetscFVCreate()
2161: */
2162: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2163: {
2164: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2165: const PetscBool useSVD = PETSC_TRUE;
2166: const PetscInt maxFaces = ls->maxFaces;
2167: PetscInt dim, f, d;
2168: PetscErrorCode ierr;
2171: if (numFaces > maxFaces) {
2172: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2173: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2174: }
2175: PetscFVGetSpatialDimension(fvm, &dim);
2176: for (f = 0; f < numFaces; ++f) {
2177: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2178: }
2179: /* Overwrites B with garbage, returns Binv in row-major format */
2180: if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2181: else {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2182: for (f = 0; f < numFaces; ++f) {
2183: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2184: }
2185: return(0);
2186: }
2188: /*
2189: neighborVol[f*2+0] contains the left geom
2190: neighborVol[f*2+1] contains the right geom
2191: */
2192: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2193: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2194: {
2195: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2196: void *rctx;
2197: PetscScalar *flux = fvm->fluxWork;
2198: const PetscScalar *constants;
2199: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2200: PetscErrorCode ierr;
2203: PetscDSGetTotalComponents(prob, &Nc);
2204: PetscDSGetTotalDimension(prob, &totDim);
2205: PetscDSGetFieldOffset(prob, field, &off);
2206: PetscDSGetRiemannSolver(prob, field, &riemann);
2207: PetscDSGetContext(prob, field, &rctx);
2208: PetscDSGetConstants(prob, &numConstants, &constants);
2209: PetscFVGetSpatialDimension(fvm, &dim);
2210: PetscFVGetNumComponents(fvm, &pdim);
2211: for (f = 0; f < Nf; ++f) {
2212: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2213: for (d = 0; d < pdim; ++d) {
2214: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2215: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2216: }
2217: }
2218: return(0);
2219: }
2221: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2222: {
2223: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2224: PetscInt dim, m, n, nrhs, minwork;
2225: PetscErrorCode ierr;
2229: PetscFVGetSpatialDimension(fvm, &dim);
2230: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2231: ls->maxFaces = maxFaces;
2232: m = ls->maxFaces;
2233: n = dim;
2234: nrhs = ls->maxFaces;
2235: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2236: ls->workSize = 5*minwork; /* We can afford to be extra generous */
2237: PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2238: return(0);
2239: }
2241: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2242: {
2244: fvm->ops->setfromoptions = NULL;
2245: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2246: fvm->ops->view = PetscFVView_LeastSquares;
2247: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2248: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2249: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2250: return(0);
2251: }
2253: /*MC
2254: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2256: Level: intermediate
2258: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2259: M*/
2261: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2262: {
2263: PetscFV_LeastSquares *ls;
2264: PetscErrorCode ierr;
2268: PetscNewLog(fvm, &ls);
2269: fvm->data = ls;
2271: ls->maxFaces = -1;
2272: ls->workSize = -1;
2273: ls->B = NULL;
2274: ls->Binv = NULL;
2275: ls->tau = NULL;
2276: ls->work = NULL;
2278: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2279: PetscFVInitialize_LeastSquares(fvm);
2280: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2281: return(0);
2282: }
2284: /*@
2285: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2287: Not collective
2289: Input parameters:
2290: + fvm - The PetscFV object
2291: - maxFaces - The maximum number of cell faces
2293: Level: intermediate
2295: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2296: @*/
2297: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2298: {
2303: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2304: return(0);
2305: }