Actual source code: dtfv.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/petscfvimpl.h> /*I "petscfv.h" I*/
2: #include <petsc/private/dmpleximpl.h> /* For CellRefiner */
3: #include <petscds.h>
5: PetscClassId PETSCLIMITER_CLASSID = 0;
7: PetscFunctionList PetscLimiterList = NULL;
8: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
10: PetscBool Limitercite = PETSC_FALSE;
11: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
12: " title = {Analysis of slope limiters on irregular grids},\n"
13: " journal = {AIAA paper},\n"
14: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
15: " volume = {490},\n"
16: " year = {2005}\n}\n";
20: /*@C
21: PetscLimiterRegister - Adds a new PetscLimiter implementation
23: Not Collective
25: Input Parameters:
26: + name - The name of a new user-defined creation routine
27: - create_func - The creation routine itself
29: Notes:
30: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
32: Sample usage:
33: .vb
34: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
35: .ve
37: Then, your PetscLimiter type can be chosen with the procedural interface via
38: .vb
39: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
40: PetscLimiterSetType(PetscLimiter, "my_lim");
41: .ve
42: or at runtime via the option
43: .vb
44: -petsclimiter_type my_lim
45: .ve
47: Level: advanced
49: .keywords: PetscLimiter, register
50: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
52: @*/
53: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
54: {
58: PetscFunctionListAdd(&PetscLimiterList, sname, function);
59: return(0);
60: }
64: /*@C
65: PetscLimiterSetType - Builds a particular PetscLimiter
67: Collective on PetscLimiter
69: Input Parameters:
70: + lim - The PetscLimiter object
71: - name - The kind of limiter
73: Options Database Key:
74: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
76: Level: intermediate
78: .keywords: PetscLimiter, set, type
79: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
80: @*/
81: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
82: {
83: PetscErrorCode (*r)(PetscLimiter);
84: PetscBool match;
89: PetscObjectTypeCompare((PetscObject) lim, name, &match);
90: if (match) return(0);
92: PetscLimiterRegisterAll();
93: PetscFunctionListFind(PetscLimiterList, name, &r);
94: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
96: if (lim->ops->destroy) {
97: (*lim->ops->destroy)(lim);
98: lim->ops->destroy = NULL;
99: }
100: (*r)(lim);
101: PetscObjectChangeTypeName((PetscObject) lim, name);
102: return(0);
103: }
107: /*@C
108: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
110: Not Collective
112: Input Parameter:
113: . lim - The PetscLimiter
115: Output Parameter:
116: . name - The PetscLimiter type name
118: Level: intermediate
120: .keywords: PetscLimiter, get, type, name
121: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
122: @*/
123: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
124: {
130: PetscLimiterRegisterAll();
131: *name = ((PetscObject) lim)->type_name;
132: return(0);
133: }
137: /*@C
138: PetscLimiterView - Views a PetscLimiter
140: Collective on PetscLimiter
142: Input Parameter:
143: + lim - the PetscLimiter object to view
144: - v - the viewer
146: Level: developer
148: .seealso: PetscLimiterDestroy()
149: @*/
150: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
151: {
156: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
157: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
158: return(0);
159: }
163: /*@
164: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
166: Collective on PetscLimiter
168: Input Parameter:
169: . lim - the PetscLimiter object to set options for
171: Level: developer
173: .seealso: PetscLimiterView()
174: @*/
175: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
176: {
177: const char *defaultType;
178: char name[256];
179: PetscBool flg;
184: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
185: else defaultType = ((PetscObject) lim)->type_name;
186: PetscLimiterRegisterAll();
188: PetscObjectOptionsBegin((PetscObject) lim);
189: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
190: if (flg) {
191: PetscLimiterSetType(lim, name);
192: } else if (!((PetscObject) lim)->type_name) {
193: PetscLimiterSetType(lim, defaultType);
194: }
195: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
196: /* process any options handlers added with PetscObjectAddOptionsHandler() */
197: PetscObjectProcessOptionsHandlers((PetscObject) lim);
198: PetscOptionsEnd();
199: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
200: return(0);
201: }
205: /*@C
206: PetscLimiterSetUp - Construct data structures for the PetscLimiter
208: Collective on PetscLimiter
210: Input Parameter:
211: . lim - the PetscLimiter object to setup
213: Level: developer
215: .seealso: PetscLimiterView(), PetscLimiterDestroy()
216: @*/
217: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
218: {
223: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
224: return(0);
225: }
229: /*@
230: PetscLimiterDestroy - Destroys a PetscLimiter object
232: Collective on PetscLimiter
234: Input Parameter:
235: . lim - the PetscLimiter object to destroy
237: Level: developer
239: .seealso: PetscLimiterView()
240: @*/
241: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
242: {
246: if (!*lim) return(0);
249: if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
250: ((PetscObject) (*lim))->refct = 0;
252: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
253: PetscHeaderDestroy(lim);
254: return(0);
255: }
259: /*@
260: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
262: Collective on MPI_Comm
264: Input Parameter:
265: . comm - The communicator for the PetscLimiter object
267: Output Parameter:
268: . lim - The PetscLimiter object
270: Level: beginner
272: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
273: @*/
274: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
275: {
276: PetscLimiter l;
281: PetscCitationsRegister(LimiterCitation,&Limitercite);
282: *lim = NULL;
283: PetscFVInitializePackage();
285: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
287: *lim = l;
288: return(0);
289: }
293: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
294: *
295: * The classical flux-limited formulation is psi(r) where
296: *
297: * r = (u[0] - u[-1]) / (u[1] - u[0])
298: *
299: * The second order TVD region is bounded by
300: *
301: * psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
302: *
303: * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
304: * phi(r)(r+1)/2 in which the reconstructed interface values are
305: *
306: * u(v) = u[0] + phi(r) (grad u)[0] v
307: *
308: * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
309: *
310: * phi_minmod(r) = 2 min(1/(1+r),r/(1+r)) phi_superbee(r) = 2 min(2/(1+r), 2r/(1+r), max(1,r)/(1+r))
311: *
312: * For a nicer symmetric formulation, rewrite in terms of
313: *
314: * f = (u[0] - u[-1]) / (u[1] - u[-1])
315: *
316: * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
317: *
318: * phi(r) = phi(1/r)
319: *
320: * becomes
321: *
322: * w(f) = w(1-f).
323: *
324: * The limiters below implement this final form w(f). The reference methods are
325: *
326: * w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
327: * */
328: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
329: {
335: (*lim->ops->limit)(lim, flim, phi);
336: return(0);
337: }
341: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
342: {
343: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
344: PetscErrorCode ierr;
347: PetscFree(l);
348: return(0);
349: }
353: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
354: {
355: PetscViewerFormat format;
356: PetscErrorCode ierr;
359: PetscViewerGetFormat(viewer, &format);
360: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
361: return(0);
362: }
366: PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
367: {
368: PetscBool iascii;
374: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
375: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
376: return(0);
377: }
381: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
382: {
384: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
385: return(0);
386: }
390: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
391: {
393: lim->ops->view = PetscLimiterView_Sin;
394: lim->ops->destroy = PetscLimiterDestroy_Sin;
395: lim->ops->limit = PetscLimiterLimit_Sin;
396: return(0);
397: }
399: /*MC
400: PETSCLIMITERSIN = "sin" - A PetscLimiter object
402: Level: intermediate
404: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
405: M*/
409: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
410: {
411: PetscLimiter_Sin *l;
412: PetscErrorCode ierr;
416: PetscNewLog(lim, &l);
417: lim->data = l;
419: PetscLimiterInitialize_Sin(lim);
420: return(0);
421: }
425: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
426: {
427: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
428: PetscErrorCode ierr;
431: PetscFree(l);
432: return(0);
433: }
437: PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
438: {
439: PetscViewerFormat format;
440: PetscErrorCode ierr;
443: PetscViewerGetFormat(viewer, &format);
444: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
445: return(0);
446: }
450: PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
451: {
452: PetscBool iascii;
458: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
459: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
460: return(0);
461: }
465: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
466: {
468: *phi = 0.0;
469: return(0);
470: }
474: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
475: {
477: lim->ops->view = PetscLimiterView_Zero;
478: lim->ops->destroy = PetscLimiterDestroy_Zero;
479: lim->ops->limit = PetscLimiterLimit_Zero;
480: return(0);
481: }
483: /*MC
484: PETSCLIMITERZERO = "zero" - A PetscLimiter object
486: Level: intermediate
488: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
489: M*/
493: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
494: {
495: PetscLimiter_Zero *l;
496: PetscErrorCode ierr;
500: PetscNewLog(lim, &l);
501: lim->data = l;
503: PetscLimiterInitialize_Zero(lim);
504: return(0);
505: }
509: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
510: {
511: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
512: PetscErrorCode ierr;
515: PetscFree(l);
516: return(0);
517: }
521: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
522: {
523: PetscViewerFormat format;
524: PetscErrorCode ierr;
527: PetscViewerGetFormat(viewer, &format);
528: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
529: return(0);
530: }
534: PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
535: {
536: PetscBool iascii;
542: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
543: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
544: return(0);
545: }
549: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
550: {
552: *phi = 1.0;
553: return(0);
554: }
558: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
559: {
561: lim->ops->view = PetscLimiterView_None;
562: lim->ops->destroy = PetscLimiterDestroy_None;
563: lim->ops->limit = PetscLimiterLimit_None;
564: return(0);
565: }
567: /*MC
568: PETSCLIMITERNONE = "none" - A PetscLimiter object
570: Level: intermediate
572: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
573: M*/
577: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
578: {
579: PetscLimiter_None *l;
580: PetscErrorCode ierr;
584: PetscNewLog(lim, &l);
585: lim->data = l;
587: PetscLimiterInitialize_None(lim);
588: return(0);
589: }
593: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
594: {
595: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
596: PetscErrorCode ierr;
599: PetscFree(l);
600: return(0);
601: }
605: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
606: {
607: PetscViewerFormat format;
608: PetscErrorCode ierr;
611: PetscViewerGetFormat(viewer, &format);
612: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
613: return(0);
614: }
618: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
619: {
620: PetscBool iascii;
626: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
627: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
628: return(0);
629: }
633: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
634: {
636: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
637: return(0);
638: }
642: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
643: {
645: lim->ops->view = PetscLimiterView_Minmod;
646: lim->ops->destroy = PetscLimiterDestroy_Minmod;
647: lim->ops->limit = PetscLimiterLimit_Minmod;
648: return(0);
649: }
651: /*MC
652: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
654: Level: intermediate
656: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
657: M*/
661: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
662: {
663: PetscLimiter_Minmod *l;
664: PetscErrorCode ierr;
668: PetscNewLog(lim, &l);
669: lim->data = l;
671: PetscLimiterInitialize_Minmod(lim);
672: return(0);
673: }
677: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
678: {
679: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
680: PetscErrorCode ierr;
683: PetscFree(l);
684: return(0);
685: }
689: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
690: {
691: PetscViewerFormat format;
692: PetscErrorCode ierr;
695: PetscViewerGetFormat(viewer, &format);
696: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
697: return(0);
698: }
702: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
703: {
704: PetscBool iascii;
710: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
711: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
712: return(0);
713: }
717: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
718: {
720: *phi = PetscMax(0, 4*f*(1-f));
721: return(0);
722: }
726: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
727: {
729: lim->ops->view = PetscLimiterView_VanLeer;
730: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
731: lim->ops->limit = PetscLimiterLimit_VanLeer;
732: return(0);
733: }
735: /*MC
736: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
738: Level: intermediate
740: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
741: M*/
745: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
746: {
747: PetscLimiter_VanLeer *l;
748: PetscErrorCode ierr;
752: PetscNewLog(lim, &l);
753: lim->data = l;
755: PetscLimiterInitialize_VanLeer(lim);
756: return(0);
757: }
761: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
762: {
763: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
764: PetscErrorCode ierr;
767: PetscFree(l);
768: return(0);
769: }
773: PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
774: {
775: PetscViewerFormat format;
776: PetscErrorCode ierr;
779: PetscViewerGetFormat(viewer, &format);
780: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
781: return(0);
782: }
786: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
787: {
788: PetscBool iascii;
794: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
795: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
796: return(0);
797: }
801: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
802: {
804: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
805: return(0);
806: }
810: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
811: {
813: lim->ops->view = PetscLimiterView_VanAlbada;
814: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
815: lim->ops->limit = PetscLimiterLimit_VanAlbada;
816: return(0);
817: }
819: /*MC
820: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
822: Level: intermediate
824: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
825: M*/
829: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
830: {
831: PetscLimiter_VanAlbada *l;
832: PetscErrorCode ierr;
836: PetscNewLog(lim, &l);
837: lim->data = l;
839: PetscLimiterInitialize_VanAlbada(lim);
840: return(0);
841: }
845: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
846: {
847: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
848: PetscErrorCode ierr;
851: PetscFree(l);
852: return(0);
853: }
857: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
858: {
859: PetscViewerFormat format;
860: PetscErrorCode ierr;
863: PetscViewerGetFormat(viewer, &format);
864: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
865: return(0);
866: }
870: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
871: {
872: PetscBool iascii;
878: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
879: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
880: return(0);
881: }
885: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
886: {
888: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
889: return(0);
890: }
894: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
895: {
897: lim->ops->view = PetscLimiterView_Superbee;
898: lim->ops->destroy = PetscLimiterDestroy_Superbee;
899: lim->ops->limit = PetscLimiterLimit_Superbee;
900: return(0);
901: }
903: /*MC
904: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
906: Level: intermediate
908: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
909: M*/
913: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
914: {
915: PetscLimiter_Superbee *l;
916: PetscErrorCode ierr;
920: PetscNewLog(lim, &l);
921: lim->data = l;
923: PetscLimiterInitialize_Superbee(lim);
924: return(0);
925: }
929: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
930: {
931: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
932: PetscErrorCode ierr;
935: PetscFree(l);
936: return(0);
937: }
941: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
942: {
943: PetscViewerFormat format;
944: PetscErrorCode ierr;
947: PetscViewerGetFormat(viewer, &format);
948: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
949: return(0);
950: }
954: PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
955: {
956: PetscBool iascii;
962: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
963: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
964: return(0);
965: }
969: /* aka Barth-Jespersen */
970: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
971: {
973: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
974: return(0);
975: }
979: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
980: {
982: lim->ops->view = PetscLimiterView_MC;
983: lim->ops->destroy = PetscLimiterDestroy_MC;
984: lim->ops->limit = PetscLimiterLimit_MC;
985: return(0);
986: }
988: /*MC
989: PETSCLIMITERMC = "mc" - A PetscLimiter object
991: Level: intermediate
993: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
994: M*/
998: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
999: {
1000: PetscLimiter_MC *l;
1001: PetscErrorCode ierr;
1005: PetscNewLog(lim, &l);
1006: lim->data = l;
1008: PetscLimiterInitialize_MC(lim);
1009: return(0);
1010: }
1012: PetscClassId PETSCFV_CLASSID = 0;
1014: PetscFunctionList PetscFVList = NULL;
1015: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
1019: /*@C
1020: PetscFVRegister - Adds a new PetscFV implementation
1022: Not Collective
1024: Input Parameters:
1025: + name - The name of a new user-defined creation routine
1026: - create_func - The creation routine itself
1028: Notes:
1029: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
1031: Sample usage:
1032: .vb
1033: PetscFVRegister("my_fv", MyPetscFVCreate);
1034: .ve
1036: Then, your PetscFV type can be chosen with the procedural interface via
1037: .vb
1038: PetscFVCreate(MPI_Comm, PetscFV *);
1039: PetscFVSetType(PetscFV, "my_fv");
1040: .ve
1041: or at runtime via the option
1042: .vb
1043: -petscfv_type my_fv
1044: .ve
1046: Level: advanced
1048: .keywords: PetscFV, register
1049: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
1051: @*/
1052: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
1053: {
1057: PetscFunctionListAdd(&PetscFVList, sname, function);
1058: return(0);
1059: }
1063: /*@C
1064: PetscFVSetType - Builds a particular PetscFV
1066: Collective on PetscFV
1068: Input Parameters:
1069: + fvm - The PetscFV object
1070: - name - The kind of FVM space
1072: Options Database Key:
1073: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
1075: Level: intermediate
1077: .keywords: PetscFV, set, type
1078: .seealso: PetscFVGetType(), PetscFVCreate()
1079: @*/
1080: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
1081: {
1082: PetscErrorCode (*r)(PetscFV);
1083: PetscBool match;
1088: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1089: if (match) return(0);
1091: PetscFVRegisterAll();
1092: PetscFunctionListFind(PetscFVList, name, &r);
1093: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1095: if (fvm->ops->destroy) {
1096: (*fvm->ops->destroy)(fvm);
1097: fvm->ops->destroy = NULL;
1098: }
1099: (*r)(fvm);
1100: PetscObjectChangeTypeName((PetscObject) fvm, name);
1101: return(0);
1102: }
1106: /*@C
1107: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1109: Not Collective
1111: Input Parameter:
1112: . fvm - The PetscFV
1114: Output Parameter:
1115: . name - The PetscFV type name
1117: Level: intermediate
1119: .keywords: PetscFV, get, type, name
1120: .seealso: PetscFVSetType(), PetscFVCreate()
1121: @*/
1122: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1123: {
1129: PetscFVRegisterAll();
1130: *name = ((PetscObject) fvm)->type_name;
1131: return(0);
1132: }
1136: /*@C
1137: PetscFVView - Views a PetscFV
1139: Collective on PetscFV
1141: Input Parameter:
1142: + fvm - the PetscFV object to view
1143: - v - the viewer
1145: Level: developer
1147: .seealso: PetscFVDestroy()
1148: @*/
1149: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1150: {
1155: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1156: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1157: return(0);
1158: }
1162: /*@
1163: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1165: Collective on PetscFV
1167: Input Parameter:
1168: . fvm - the PetscFV object to set options for
1170: Level: developer
1172: .seealso: PetscFVView()
1173: @*/
1174: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1175: {
1176: const char *defaultType;
1177: char name[256];
1178: PetscBool flg;
1183: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1184: else defaultType = ((PetscObject) fvm)->type_name;
1185: PetscFVRegisterAll();
1187: PetscObjectOptionsBegin((PetscObject) fvm);
1188: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1189: if (flg) {
1190: PetscFVSetType(fvm, name);
1191: } else if (!((PetscObject) fvm)->type_name) {
1192: PetscFVSetType(fvm, defaultType);
1193: }
1194: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1195: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1196: PetscObjectProcessOptionsHandlers((PetscObject) fvm);
1197: PetscLimiterSetFromOptions(fvm->limiter);
1198: PetscOptionsEnd();
1199: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1200: return(0);
1201: }
1205: /*@
1206: PetscFVSetUp - Construct data structures for the PetscFV
1208: Collective on PetscFV
1210: Input Parameter:
1211: . fvm - the PetscFV object to setup
1213: Level: developer
1215: .seealso: PetscFVView(), PetscFVDestroy()
1216: @*/
1217: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1218: {
1223: PetscLimiterSetUp(fvm->limiter);
1224: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1225: return(0);
1226: }
1230: /*@
1231: PetscFVDestroy - Destroys a PetscFV object
1233: Collective on PetscFV
1235: Input Parameter:
1236: . fvm - the PetscFV object to destroy
1238: Level: developer
1240: .seealso: PetscFVView()
1241: @*/
1242: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1243: {
1247: if (!*fvm) return(0);
1250: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1251: ((PetscObject) (*fvm))->refct = 0;
1253: PetscLimiterDestroy(&(*fvm)->limiter);
1254: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1255: PetscFree((*fvm)->fluxWork);
1256: PetscQuadratureDestroy(&(*fvm)->quadrature);
1257: PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);
1259: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1260: PetscHeaderDestroy(fvm);
1261: return(0);
1262: }
1266: /*@
1267: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1269: Collective on MPI_Comm
1271: Input Parameter:
1272: . comm - The communicator for the PetscFV object
1274: Output Parameter:
1275: . fvm - The PetscFV object
1277: Level: beginner
1279: .seealso: PetscFVSetType(), PETSCFVUPWIND
1280: @*/
1281: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1282: {
1283: PetscFV f;
1288: *fvm = NULL;
1289: PetscFVInitializePackage();
1291: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1292: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1294: PetscLimiterCreate(comm, &f->limiter);
1295: f->numComponents = 1;
1296: f->dim = 0;
1297: f->computeGradients = PETSC_FALSE;
1298: f->fluxWork = NULL;
1300: *fvm = f;
1301: return(0);
1302: }
1306: /*@
1307: PetscFVSetLimiter - Set the limiter object
1309: Logically collective on PetscFV
1311: Input Parameters:
1312: + fvm - the PetscFV object
1313: - lim - The PetscLimiter
1315: Level: developer
1317: .seealso: PetscFVGetLimiter()
1318: @*/
1319: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1320: {
1326: PetscLimiterDestroy(&fvm->limiter);
1327: PetscObjectReference((PetscObject) lim);
1328: fvm->limiter = lim;
1329: return(0);
1330: }
1334: /*@
1335: PetscFVGetLimiter - Get the limiter object
1337: Not collective
1339: Input Parameter:
1340: . fvm - the PetscFV object
1342: Output Parameter:
1343: . lim - The PetscLimiter
1345: Level: developer
1347: .seealso: PetscFVSetLimiter()
1348: @*/
1349: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1350: {
1354: *lim = fvm->limiter;
1355: return(0);
1356: }
1360: /*@
1361: PetscFVSetNumComponents - Set the number of field components
1363: Logically collective on PetscFV
1365: Input Parameters:
1366: + fvm - the PetscFV object
1367: - comp - The number of components
1369: Level: developer
1371: .seealso: PetscFVGetNumComponents()
1372: @*/
1373: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1374: {
1379: fvm->numComponents = comp;
1380: PetscFree(fvm->fluxWork);
1381: PetscMalloc1(comp, &fvm->fluxWork);
1382: return(0);
1383: }
1387: /*@
1388: PetscFVGetNumComponents - Get the number of field components
1390: Not collective
1392: Input Parameter:
1393: . fvm - the PetscFV object
1395: Output Parameter:
1396: , comp - The number of components
1398: Level: developer
1400: .seealso: PetscFVSetNumComponents()
1401: @*/
1402: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1403: {
1407: *comp = fvm->numComponents;
1408: return(0);
1409: }
1413: /*@
1414: PetscFVSetSpatialDimension - Set the spatial dimension
1416: Logically collective on PetscFV
1418: Input Parameters:
1419: + fvm - the PetscFV object
1420: - dim - The spatial dimension
1422: Level: developer
1424: .seealso: PetscFVGetSpatialDimension()
1425: @*/
1426: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1427: {
1430: fvm->dim = dim;
1431: return(0);
1432: }
1436: /*@
1437: PetscFVGetSpatialDimension - Get the spatial dimension
1439: Logically collective on PetscFV
1441: Input Parameter:
1442: . fvm - the PetscFV object
1444: Output Parameter:
1445: . dim - The spatial dimension
1447: Level: developer
1449: .seealso: PetscFVSetSpatialDimension()
1450: @*/
1451: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1452: {
1456: *dim = fvm->dim;
1457: return(0);
1458: }
1462: /*@
1463: PetscFVSetComputeGradients - Toggle computation of cell gradients
1465: Logically collective on PetscFV
1467: Input Parameters:
1468: + fvm - the PetscFV object
1469: - computeGradients - Flag to compute cell gradients
1471: Level: developer
1473: .seealso: PetscFVGetComputeGradients()
1474: @*/
1475: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1476: {
1479: fvm->computeGradients = computeGradients;
1480: return(0);
1481: }
1485: /*@
1486: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1488: Not collective
1490: Input Parameter:
1491: . fvm - the PetscFV object
1493: Output Parameter:
1494: . computeGradients - Flag to compute cell gradients
1496: Level: developer
1498: .seealso: PetscFVSetComputeGradients()
1499: @*/
1500: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1501: {
1505: *computeGradients = fvm->computeGradients;
1506: return(0);
1507: }
1511: /*@
1512: PetscFVSetQuadrature - Set the quadrature object
1514: Logically collective on PetscFV
1516: Input Parameters:
1517: + fvm - the PetscFV object
1518: - q - The PetscQuadrature
1520: Level: developer
1522: .seealso: PetscFVGetQuadrature()
1523: @*/
1524: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1525: {
1530: PetscQuadratureDestroy(&fvm->quadrature);
1531: PetscObjectReference((PetscObject) q);
1532: fvm->quadrature = q;
1533: return(0);
1534: }
1538: /*@
1539: PetscFVGetQuadrature - Get the quadrature object
1541: Not collective
1543: Input Parameter:
1544: . fvm - the PetscFV object
1546: Output Parameter:
1547: . lim - The PetscQuadrature
1549: Level: developer
1551: .seealso: PetscFVSetQuadrature()
1552: @*/
1553: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1554: {
1558: if (!fvm->quadrature) {
1559: /* Create default 1-point quadrature */
1560: PetscReal *points, *weights;
1563: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1564: PetscCalloc1(fvm->dim, &points);
1565: PetscMalloc1(1, &weights);
1566: weights[0] = 1.0;
1567: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, points, weights);
1568: }
1569: *q = fvm->quadrature;
1570: return(0);
1571: }
1575: /*@
1576: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1578: Not collective
1580: Input Parameter:
1581: . fvm - The PetscFV object
1583: Output Parameter:
1584: . sp - The PetscDualSpace object
1586: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1588: Level: developer
1590: .seealso: PetscFVCreate()
1591: @*/
1592: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1593: {
1597: if (!fvm->dualSpace) {
1598: DM K;
1599: PetscQuadrature q;
1600: PetscInt dim;
1601: PetscErrorCode ierr;
1603: PetscFVGetSpatialDimension(fvm, &dim);
1604: PetscFVGetQuadrature(fvm, &q);
1605: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1606: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1607: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, 1);
1608: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, 0, q);
1609: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1610: PetscDualSpaceSetDM(fvm->dualSpace, K);
1611: DMDestroy(&K);
1612: }
1613: *sp = fvm->dualSpace;
1614: return(0);
1615: }
1619: /*@
1620: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1622: Not collective
1624: Input Parameters:
1625: + fvm - The PetscFV object
1626: - sp - The PetscDualSpace object
1628: Level: intermediate
1630: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1632: .seealso: PetscFVCreate()
1633: @*/
1634: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1635: {
1641: PetscDualSpaceDestroy(&fvm->dualSpace);
1642: fvm->dualSpace = sp;
1643: PetscObjectReference((PetscObject) fvm->dualSpace);
1644: return(0);
1645: }
1649: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1650: {
1651: PetscInt npoints;
1652: const PetscReal *points;
1653: PetscErrorCode ierr;
1660: PetscQuadratureGetData(fvm->quadrature, NULL, &npoints, &points, NULL);
1661: if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1662: if (B) *B = fvm->B;
1663: if (D) *D = fvm->D;
1664: if (H) *H = fvm->H;
1665: return(0);
1666: }
1670: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1671: {
1672: PetscInt pdim = 1; /* Dimension of approximation space P */
1673: PetscInt dim; /* Spatial dimension */
1674: PetscInt comp; /* Field components */
1675: PetscInt p, d, c, e;
1676: PetscErrorCode ierr;
1684: PetscFVGetSpatialDimension(fvm, &dim);
1685: PetscFVGetNumComponents(fvm, &comp);
1686: if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1687: if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1688: if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1689: 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;}
1690: 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;}
1691: 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;}
1692: return(0);
1693: }
1697: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1698: {
1703: if (B && *B) {PetscFree(*B);}
1704: if (D && *D) {PetscFree(*D);}
1705: if (H && *H) {PetscFree(*H);}
1706: return(0);
1707: }
1711: /*@C
1712: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1714: Input Parameters:
1715: + fvm - The PetscFV object
1716: . numFaces - The number of cell faces which are not constrained
1717: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1719: Level: developer
1721: .seealso: PetscFVCreate()
1722: @*/
1723: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1724: {
1729: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1730: return(0);
1731: }
1735: /*C
1736: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1738: Not collective
1740: Input Parameters:
1741: + fvm - The PetscFV object for the field being integrated
1742: . prob - The PetscDS specifing the discretizations and continuum functions
1743: . field - The field being integrated
1744: . Nf - The number of faces in the chunk
1745: . fgeom - The face geometry for each face in the chunk
1746: . neighborVol - The volume for each pair of cells in the chunk
1747: . uL - The state from the cell on the left
1748: - uR - The state from the cell on the right
1750: Output Parameter
1751: + fluxL - the left fluxes for each face
1752: - fluxR - the right fluxes for each face
1753: */
1754: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1755: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1756: {
1761: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1762: return(0);
1763: }
1767: /*@
1768: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1769: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1770: sparsity). It is also used to create an interpolation between regularly refined meshes.
1772: Input Parameter:
1773: . fv - The initial PetscFV
1775: Output Parameter:
1776: . fvRef - The refined PetscFV
1778: Level: developer
1780: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1781: @*/
1782: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1783: {
1784: PetscDualSpace Q, Qref;
1785: DM K, Kref;
1786: PetscQuadrature q, qref;
1787: CellRefiner cellRefiner;
1788: PetscReal *v0;
1789: PetscReal *jac, *invjac;
1790: PetscInt numComp, numSubelements, s;
1791: PetscErrorCode ierr;
1794: PetscFVGetDualSpace(fv, &Q);
1795: PetscFVGetQuadrature(fv, &q);
1796: PetscDualSpaceGetDM(Q, &K);
1797: /* Create dual space */
1798: PetscDualSpaceDuplicate(Q, &Qref);
1799: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1800: PetscDualSpaceSetDM(Qref, Kref);
1801: DMDestroy(&Kref);
1802: PetscDualSpaceSetUp(Qref);
1803: /* Create volume */
1804: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1805: PetscFVSetDualSpace(*fvRef, Qref);
1806: PetscFVGetNumComponents(fv, &numComp);
1807: PetscFVSetNumComponents(*fvRef, numComp);
1808: PetscFVSetUp(*fvRef);
1809: /* Create quadrature */
1810: DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1811: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1812: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1813: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1814: for (s = 0; s < numSubelements; ++s) {
1815: PetscQuadrature qs;
1816: const PetscReal *points, *weights;
1817: PetscReal *p, *w;
1818: PetscInt dim, npoints, np;
1820: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1821: PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
1822: np = npoints/numSubelements;
1823: PetscMalloc1(np*dim,&p);
1824: PetscMalloc1(np,&w);
1825: PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1826: PetscMemcpy(w, &weights[s*np], np * sizeof(PetscReal));
1827: PetscQuadratureSetData(qs, dim, np, p, w);
1828: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1829: PetscQuadratureDestroy(&qs);
1830: }
1831: CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1832: PetscFVSetQuadrature(*fvRef, qref);
1833: PetscQuadratureDestroy(&qref);
1834: PetscDualSpaceDestroy(&Qref);
1835: return(0);
1836: }
1840: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1841: {
1842: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1846: PetscFree(b);
1847: return(0);
1848: }
1852: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1853: {
1854: PetscInt Nc;
1855: PetscViewerFormat format;
1856: PetscErrorCode ierr;
1859: PetscFVGetNumComponents(fv, &Nc);
1860: PetscViewerGetFormat(viewer, &format);
1861: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1862: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1863: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1864: } else {
1865: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1866: }
1867: return(0);
1868: }
1872: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1873: {
1874: PetscBool iascii;
1880: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1881: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1882: return(0);
1883: }
1887: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1888: {
1890: return(0);
1891: }
1895: /*
1896: neighborVol[f*2+0] contains the left geom
1897: neighborVol[f*2+1] contains the right geom
1898: */
1899: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1900: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1901: {
1902: void (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
1903: void *rctx;
1904: PetscScalar *flux = fvm->fluxWork;
1905: PetscInt dim, pdim, totDim, Nc, off, f, d;
1909: PetscDSGetTotalComponents(prob, &Nc);
1910: PetscDSGetTotalDimension(prob, &totDim);
1911: PetscDSGetFieldOffset(prob, field, &off);
1912: PetscDSGetRiemannSolver(prob, field, &riemann);
1913: PetscDSGetContext(prob, field, &rctx);
1914: PetscFVGetSpatialDimension(fvm, &dim);
1915: PetscFVGetNumComponents(fvm, &pdim);
1916: for (f = 0; f < Nf; ++f) {
1917: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
1918: for (d = 0; d < pdim; ++d) {
1919: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1920: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1921: }
1922: }
1923: return(0);
1924: }
1928: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1929: {
1931: fvm->ops->setfromoptions = NULL;
1932: fvm->ops->setup = PetscFVSetUp_Upwind;
1933: fvm->ops->view = PetscFVView_Upwind;
1934: fvm->ops->destroy = PetscFVDestroy_Upwind;
1935: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1936: return(0);
1937: }
1939: /*MC
1940: PETSCFVUPWIND = "upwind" - A PetscFV object
1942: Level: intermediate
1944: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1945: M*/
1949: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1950: {
1951: PetscFV_Upwind *b;
1952: PetscErrorCode ierr;
1956: PetscNewLog(fvm,&b);
1957: fvm->data = b;
1959: PetscFVInitialize_Upwind(fvm);
1960: return(0);
1961: }
1963: #include <petscblaslapack.h>
1967: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1968: {
1969: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1970: PetscErrorCode ierr;
1973: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1974: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1975: PetscFree(ls);
1976: return(0);
1977: }
1981: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1982: {
1983: PetscInt Nc;
1984: PetscViewerFormat format;
1985: PetscErrorCode ierr;
1988: PetscFVGetNumComponents(fv, &Nc);
1989: PetscViewerGetFormat(viewer, &format);
1990: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1991: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1992: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1993: } else {
1994: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1995: }
1996: return(0);
1997: }
2001: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2002: {
2003: PetscBool iascii;
2009: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2010: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2011: return(0);
2012: }
2016: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2017: {
2019: return(0);
2020: }
2024: /* Overwrites A. Can only handle full-rank problems with m>=n */
2025: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2026: {
2027: PetscBool debug = PETSC_FALSE;
2029: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2030: PetscScalar *R,*Q,*Aback,Alpha;
2033: if (debug) {
2034: PetscMalloc1(m*n,&Aback);
2035: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2036: }
2038: PetscBLASIntCast(m,&M);
2039: PetscBLASIntCast(n,&N);
2040: PetscBLASIntCast(mstride,&lda);
2041: PetscBLASIntCast(worksize,&ldwork);
2042: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2043: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2044: PetscFPTrapPop();
2045: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2046: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2048: /* Extract an explicit representation of Q */
2049: Q = Ainv;
2050: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
2051: K = N; /* full rank */
2052: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
2053: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2055: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2056: Alpha = 1.0;
2057: ldb = lda;
2058: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2059: /* Ainv is Q, overwritten with inverse */
2061: if (debug) { /* Check that pseudo-inverse worked */
2062: PetscScalar Beta = 0.0;
2063: PetscBLASInt ldc;
2064: K = N;
2065: ldc = N;
2066: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2067: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2068: PetscFree(Aback);
2069: }
2070: return(0);
2071: }
2075: /* Overwrites A. Can handle degenerate problems and m<n. */
2076: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2077: {
2078: PetscBool debug = PETSC_FALSE;
2079: PetscScalar *Brhs, *Aback;
2080: PetscScalar *tmpwork;
2081: PetscReal rcond;
2082: PetscInt i, j, maxmn;
2083: PetscBLASInt M, N, lda, ldb, ldwork;
2084: #if !defined(PETSC_USE_COMPLEX)
2085: PetscBLASInt nrhs, irank, info;
2086: #endif
2090: if (debug) {
2091: PetscMalloc1(m*n,&Aback);
2092: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2093: }
2095: /* initialize to identity */
2096: tmpwork = Ainv;
2097: Brhs = work;
2098: maxmn = PetscMax(m,n);
2099: for (j=0; j<maxmn; j++) {
2100: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2101: }
2103: PetscBLASIntCast(m,&M);
2104: PetscBLASIntCast(n,&N);
2105: PetscBLASIntCast(mstride,&lda);
2106: PetscBLASIntCast(maxmn,&ldb);
2107: PetscBLASIntCast(worksize,&ldwork);
2108: rcond = -1;
2109: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2110: #if defined(PETSC_USE_COMPLEX)
2111: if (tmpwork && rcond) rcond = 0.0; /* Get rid of compiler warning */
2112: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "I don't think this makes sense for complex numbers");
2113: #else
2114: nrhs = M;
2115: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2116: PetscFPTrapPop();
2117: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2118: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2119: 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");
2121: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2122: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2123: for (i=0; i<n; i++) {
2124: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2125: }
2126: return(0);
2127: #endif
2128: }
2130: #if 0
2133: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2134: {
2135: PetscReal grad[2] = {0, 0};
2136: const PetscInt *faces;
2137: PetscInt numFaces, f;
2138: PetscErrorCode ierr;
2141: DMPlexGetConeSize(dm, cell, &numFaces);
2142: DMPlexGetCone(dm, cell, &faces);
2143: for (f = 0; f < numFaces; ++f) {
2144: const PetscInt *fcells;
2145: const CellGeom *cg1;
2146: const FaceGeom *fg;
2148: DMPlexGetSupport(dm, faces[f], &fcells);
2149: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2150: for (i = 0; i < 2; ++i) {
2151: PetscScalar du;
2153: if (fcells[i] == c) continue;
2154: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2155: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2156: grad[0] += fg->grad[!i][0] * du;
2157: grad[1] += fg->grad[!i][1] * du;
2158: }
2159: }
2160: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2161: return(0);
2162: }
2163: #endif
2167: /*
2168: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2170: Input Parameters:
2171: + fvm - The PetscFV object
2172: . numFaces - The number of cell faces which are not constrained
2173: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2175: Level: developer
2177: .seealso: PetscFVCreate()
2178: */
2179: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2180: {
2181: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2182: const PetscBool useSVD = PETSC_TRUE;
2183: const PetscInt maxFaces = ls->maxFaces;
2184: PetscInt dim, f, d;
2185: PetscErrorCode ierr;
2188: if (numFaces > maxFaces) {
2189: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2190: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2191: }
2192: PetscFVGetSpatialDimension(fvm, &dim);
2193: for (f = 0; f < numFaces; ++f) {
2194: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2195: }
2196: /* Overwrites B with garbage, returns Binv in row-major format */
2197: if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2198: else {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2199: for (f = 0; f < numFaces; ++f) {
2200: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2201: }
2202: return(0);
2203: }
2207: /*
2208: neighborVol[f*2+0] contains the left geom
2209: neighborVol[f*2+1] contains the right geom
2210: */
2211: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2212: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2213: {
2214: void (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
2215: void *rctx;
2216: PetscScalar *flux = fvm->fluxWork;
2217: PetscInt dim, pdim, Nc, totDim, off, f, d;
2221: PetscDSGetTotalComponents(prob, &Nc);
2222: PetscDSGetTotalDimension(prob, &totDim);
2223: PetscDSGetFieldOffset(prob, field, &off);
2224: PetscDSGetRiemannSolver(prob, field, &riemann);
2225: PetscDSGetContext(prob, field, &rctx);
2226: PetscFVGetSpatialDimension(fvm, &dim);
2227: PetscFVGetNumComponents(fvm, &pdim);
2228: for (f = 0; f < Nf; ++f) {
2229: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
2230: for (d = 0; d < pdim; ++d) {
2231: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2232: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2233: }
2234: }
2235: return(0);
2236: }
2240: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2241: {
2242: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2243: PetscInt dim, m, n, nrhs, minwork;
2244: PetscErrorCode ierr;
2248: PetscFVGetSpatialDimension(fvm, &dim);
2249: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2250: ls->maxFaces = maxFaces;
2251: m = ls->maxFaces;
2252: n = dim;
2253: nrhs = ls->maxFaces;
2254: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2255: ls->workSize = 5*minwork; /* We can afford to be extra generous */
2256: PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2257: return(0);
2258: }
2262: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2263: {
2265: fvm->ops->setfromoptions = NULL;
2266: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2267: fvm->ops->view = PetscFVView_LeastSquares;
2268: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2269: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2270: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2271: return(0);
2272: }
2274: /*MC
2275: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2277: Level: intermediate
2279: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2280: M*/
2284: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2285: {
2286: PetscFV_LeastSquares *ls;
2287: PetscErrorCode ierr;
2291: PetscNewLog(fvm, &ls);
2292: fvm->data = ls;
2294: ls->maxFaces = -1;
2295: ls->workSize = -1;
2296: ls->B = NULL;
2297: ls->Binv = NULL;
2298: ls->tau = NULL;
2299: ls->work = NULL;
2301: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2302: PetscFVInitialize_LeastSquares(fvm);
2303: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2304: return(0);
2305: }
2309: /*@
2310: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2312: Not collective
2314: Input parameters:
2315: + fvm - The PetscFV object
2316: - maxFaces - The maximum number of cell faces
2318: Level: intermediate
2320: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2321: @*/
2322: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2323: {
2328: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2329: return(0);
2330: }