Actual source code: dtfv.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/petscfvimpl.h> /*I "petscfv.h" I*/
2: #include <petscdmplex.h>
4: PetscClassId PETSCLIMITER_CLASSID = 0;
6: PetscFunctionList PetscLimiterList = NULL;
7: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
9: PetscBool Limitercite = PETSC_FALSE;
10: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
11: " title = {Analysis of slope limiters on irregular grids},\n"
12: " journal = {AIAA paper},\n"
13: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
14: " volume = {490},\n"
15: " year = {2005}\n}\n";
19: /*@C
20: PetscLimiterRegister - Adds a new PetscLimiter implementation
22: Not Collective
24: Input Parameters:
25: + name - The name of a new user-defined creation routine
26: - create_func - The creation routine itself
28: Notes:
29: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
31: Sample usage:
32: .vb
33: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
34: .ve
36: Then, your PetscLimiter type can be chosen with the procedural interface via
37: .vb
38: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
39: PetscLimiterSetType(PetscLimiter, "my_lim");
40: .ve
41: or at runtime via the option
42: .vb
43: -petsclimiter_type my_lim
44: .ve
46: Level: advanced
48: .keywords: PetscLimiter, register
49: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
51: @*/
52: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
53: {
57: PetscFunctionListAdd(&PetscLimiterList, sname, function);
58: return(0);
59: }
63: /*@C
64: PetscLimiterSetType - Builds a particular PetscLimiter
66: Collective on PetscLimiter
68: Input Parameters:
69: + lim - The PetscLimiter object
70: - name - The kind of limiter
72: Options Database Key:
73: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
75: Level: intermediate
77: .keywords: PetscLimiter, set, type
78: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
79: @*/
80: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
81: {
82: PetscErrorCode (*r)(PetscLimiter);
83: PetscBool match;
88: PetscObjectTypeCompare((PetscObject) lim, name, &match);
89: if (match) return(0);
91: if (!PetscLimiterRegisterAllCalled) {PetscLimiterRegisterAll();}
92: PetscFunctionListFind(PetscLimiterList, name, &r);
93: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
95: if (lim->ops->destroy) {
96: (*lim->ops->destroy)(lim);
97: lim->ops->destroy = NULL;
98: }
99: (*r)(lim);
100: PetscObjectChangeTypeName((PetscObject) lim, name);
101: return(0);
102: }
106: /*@C
107: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
109: Not Collective
111: Input Parameter:
112: . lim - The PetscLimiter
114: Output Parameter:
115: . name - The PetscLimiter type name
117: Level: intermediate
119: .keywords: PetscLimiter, get, type, name
120: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
121: @*/
122: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
123: {
129: if (!PetscLimiterRegisterAllCalled) {PetscLimiterRegisterAll();}
130: *name = ((PetscObject) lim)->type_name;
131: return(0);
132: }
136: /*@C
137: PetscLimiterView - Views a PetscLimiter
139: Collective on PetscLimiter
141: Input Parameter:
142: + lim - the PetscLimiter object to view
143: - v - the viewer
145: Level: developer
147: .seealso: PetscLimiterDestroy()
148: @*/
149: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
150: {
155: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
156: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
157: return(0);
158: }
162: /*
163: PetscLimiterViewFromOptions - Processes command line options to determine if/how a PetscLimiter is to be viewed.
165: Collective on PetscLimiter
167: Input Parameters:
168: + lim - the PetscLimiter
169: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
170: - optionname - option to activate viewing
172: Level: intermediate
174: .keywords: PetscLimiter, view, options, database
175: .seealso: VecViewFromOptions(), MatViewFromOptions()
176: */
177: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter lim, const char prefix[], const char optionname[])
178: {
179: PetscViewer viewer;
180: PetscViewerFormat format;
181: PetscBool flg;
182: PetscErrorCode ierr;
185: if (prefix) {PetscOptionsGetViewer(PetscObjectComm((PetscObject) lim), prefix, optionname, &viewer, &format, &flg);}
186: else {PetscOptionsGetViewer(PetscObjectComm((PetscObject) lim), ((PetscObject) lim)->prefix, optionname, &viewer, &format, &flg);}
187: if (flg) {
188: PetscViewerPushFormat(viewer, format);
189: PetscLimiterView(lim, viewer);
190: PetscViewerPopFormat(viewer);
191: PetscViewerDestroy(&viewer);
192: }
193: return(0);
194: }
198: /*@
199: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
201: Collective on PetscLimiter
203: Input Parameter:
204: . lim - the PetscLimiter object to set options for
206: Level: developer
208: .seealso: PetscLimiterView()
209: @*/
210: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
211: {
212: const char *defaultType;
213: char name[256];
214: PetscBool flg;
219: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
220: else defaultType = ((PetscObject) lim)->type_name;
221: if (!PetscLimiterRegisterAllCalled) {PetscLimiterRegisterAll();}
223: PetscObjectOptionsBegin((PetscObject) lim);
224: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
225: if (flg) {
226: PetscLimiterSetType(lim, name);
227: } else if (!((PetscObject) lim)->type_name) {
228: PetscLimiterSetType(lim, defaultType);
229: }
230: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
231: /* process any options handlers added with PetscObjectAddOptionsHandler() */
232: PetscObjectProcessOptionsHandlers((PetscObject) lim);
233: PetscOptionsEnd();
234: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
235: return(0);
236: }
240: /*@C
241: PetscLimiterSetUp - Construct data structures for the PetscLimiter
243: Collective on PetscLimiter
245: Input Parameter:
246: . lim - the PetscLimiter object to setup
248: Level: developer
250: .seealso: PetscLimiterView(), PetscLimiterDestroy()
251: @*/
252: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
253: {
258: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
259: return(0);
260: }
264: /*@
265: PetscLimiterDestroy - Destroys a PetscLimiter object
267: Collective on PetscLimiter
269: Input Parameter:
270: . lim - the PetscLimiter object to destroy
272: Level: developer
274: .seealso: PetscLimiterView()
275: @*/
276: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
277: {
281: if (!*lim) return(0);
284: if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
285: ((PetscObject) (*lim))->refct = 0;
287: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
288: PetscHeaderDestroy(lim);
289: return(0);
290: }
294: /*@
295: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
297: Collective on MPI_Comm
299: Input Parameter:
300: . comm - The communicator for the PetscLimiter object
302: Output Parameter:
303: . lim - The PetscLimiter object
305: Level: beginner
307: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
308: @*/
309: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
310: {
311: PetscLimiter l;
316: PetscCitationsRegister(LimiterCitation,&Limitercite);
317: *lim = NULL;
318: PetscFVInitializePackage();
320: PetscHeaderCreate(l, _p_PetscLimiter, struct _PetscLimiterOps, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
321: PetscMemzero(l->ops, sizeof(struct _PetscLimiterOps));
323: *lim = l;
324: return(0);
325: }
329: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
330: *
331: * The classical flux-limited formulation is psi(r) where
332: *
333: * r = (u[0] - u[-1]) / (u[1] - u[0])
334: *
335: * The second order TVD region is bounded by
336: *
337: * psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
338: *
339: * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
340: * phi(r)(r+1)/2 in which the reconstructed interface values are
341: *
342: * u(v) = u[0] + phi(r) (grad u)[0] v
343: *
344: * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
345: *
346: * 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))
347: *
348: * For a nicer symmetric formulation, rewrite in terms of
349: *
350: * f = (u[0] - u[-1]) / (u[1] - u[-1])
351: *
352: * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
353: *
354: * phi(r) = phi(1/r)
355: *
356: * becomes
357: *
358: * w(f) = w(1-f).
359: *
360: * The limiters below implement this final form w(f). The reference methods are
361: *
362: * w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
363: * */
364: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
365: {
371: (*lim->ops->limit)(lim, flim, phi);
372: return(0);
373: }
377: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter fvm)
378: {
379: PetscLimiter_Sin *l = (PetscLimiter_Sin *) fvm->data;
380: PetscErrorCode ierr;
383: PetscFree(l);
384: return(0);
385: }
389: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter fv, PetscViewer viewer)
390: {
391: PetscViewerFormat format;
392: PetscErrorCode ierr;
395: PetscViewerGetFormat(viewer, &format);
396: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
397: return(0);
398: }
402: PetscErrorCode PetscLimiterView_Sin(PetscLimiter fv, PetscViewer viewer)
403: {
404: PetscBool iascii;
410: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
411: if (iascii) {PetscLimiterView_Sin_Ascii(fv, viewer);}
412: return(0);
413: }
417: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
418: {
420: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
421: return(0);
422: }
426: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter fvm)
427: {
429: fvm->ops->view = PetscLimiterView_Sin;
430: fvm->ops->destroy = PetscLimiterDestroy_Sin;
431: fvm->ops->limit = PetscLimiterLimit_Sin;
432: return(0);
433: }
435: /*MC
436: PETSCLIMITERSIN = "sin" - A PetscLimiter object
438: Level: intermediate
440: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
441: M*/
445: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
446: {
447: PetscLimiter_Sin *l;
448: PetscErrorCode ierr;
452: PetscNewLog(lim, &l);
453: lim->data = l;
455: PetscLimiterInitialize_Sin(lim);
456: return(0);
457: }
461: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter fvm)
462: {
463: PetscLimiter_Zero *l = (PetscLimiter_Zero *) fvm->data;
464: PetscErrorCode ierr;
467: PetscFree(l);
468: return(0);
469: }
473: PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter fv, PetscViewer viewer)
474: {
475: PetscViewerFormat format;
476: PetscErrorCode ierr;
479: PetscViewerGetFormat(viewer, &format);
480: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
481: return(0);
482: }
486: PetscErrorCode PetscLimiterView_Zero(PetscLimiter fv, PetscViewer viewer)
487: {
488: PetscBool iascii;
494: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
495: if (iascii) {PetscLimiterView_Zero_Ascii(fv, viewer);}
496: return(0);
497: }
501: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
502: {
504: *phi = 0.0;
505: return(0);
506: }
510: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter fvm)
511: {
513: fvm->ops->view = PetscLimiterView_Zero;
514: fvm->ops->destroy = PetscLimiterDestroy_Zero;
515: fvm->ops->limit = PetscLimiterLimit_Zero;
516: return(0);
517: }
519: /*MC
520: PETSCLIMITERZERO = "zero" - A PetscLimiter object
522: Level: intermediate
524: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
525: M*/
529: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
530: {
531: PetscLimiter_Zero *l;
532: PetscErrorCode ierr;
536: PetscNewLog(lim, &l);
537: lim->data = l;
539: PetscLimiterInitialize_Zero(lim);
540: return(0);
541: }
545: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter fvm)
546: {
547: PetscLimiter_None *l = (PetscLimiter_None *) fvm->data;
548: PetscErrorCode ierr;
551: PetscFree(l);
552: return(0);
553: }
557: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter fv, PetscViewer viewer)
558: {
559: PetscViewerFormat format;
560: PetscErrorCode ierr;
563: PetscViewerGetFormat(viewer, &format);
564: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
565: return(0);
566: }
570: PetscErrorCode PetscLimiterView_None(PetscLimiter fv, PetscViewer viewer)
571: {
572: PetscBool iascii;
578: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
579: if (iascii) {PetscLimiterView_None_Ascii(fv, viewer);}
580: return(0);
581: }
585: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
586: {
588: *phi = 1.0;
589: return(0);
590: }
594: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter fvm)
595: {
597: fvm->ops->view = PetscLimiterView_None;
598: fvm->ops->destroy = PetscLimiterDestroy_None;
599: fvm->ops->limit = PetscLimiterLimit_None;
600: return(0);
601: }
603: /*MC
604: PETSCLIMITERNONE = "none" - A PetscLimiter object
606: Level: intermediate
608: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
609: M*/
613: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
614: {
615: PetscLimiter_None *l;
616: PetscErrorCode ierr;
620: PetscNewLog(lim, &l);
621: lim->data = l;
623: PetscLimiterInitialize_None(lim);
624: return(0);
625: }
629: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter fvm)
630: {
631: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) fvm->data;
632: PetscErrorCode ierr;
635: PetscFree(l);
636: return(0);
637: }
641: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter fv, PetscViewer viewer)
642: {
643: PetscViewerFormat format;
644: PetscErrorCode ierr;
647: PetscViewerGetFormat(viewer, &format);
648: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
649: return(0);
650: }
654: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter fv, PetscViewer viewer)
655: {
656: PetscBool iascii;
662: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
663: if (iascii) {PetscLimiterView_Minmod_Ascii(fv, viewer);}
664: return(0);
665: }
669: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
670: {
672: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
673: return(0);
674: }
678: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter fvm)
679: {
681: fvm->ops->view = PetscLimiterView_Minmod;
682: fvm->ops->destroy = PetscLimiterDestroy_Minmod;
683: fvm->ops->limit = PetscLimiterLimit_Minmod;
684: return(0);
685: }
687: /*MC
688: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
690: Level: intermediate
692: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
693: M*/
697: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
698: {
699: PetscLimiter_Minmod *l;
700: PetscErrorCode ierr;
704: PetscNewLog(lim, &l);
705: lim->data = l;
707: PetscLimiterInitialize_Minmod(lim);
708: return(0);
709: }
713: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter fvm)
714: {
715: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) fvm->data;
716: PetscErrorCode ierr;
719: PetscFree(l);
720: return(0);
721: }
725: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter fv, PetscViewer viewer)
726: {
727: PetscViewerFormat format;
728: PetscErrorCode ierr;
731: PetscViewerGetFormat(viewer, &format);
732: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
733: return(0);
734: }
738: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter fv, PetscViewer viewer)
739: {
740: PetscBool iascii;
746: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
747: if (iascii) {PetscLimiterView_VanLeer_Ascii(fv, viewer);}
748: return(0);
749: }
753: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
754: {
756: *phi = PetscMax(0, 4*f*(1-f));
757: return(0);
758: }
762: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter fvm)
763: {
765: fvm->ops->view = PetscLimiterView_VanLeer;
766: fvm->ops->destroy = PetscLimiterDestroy_VanLeer;
767: fvm->ops->limit = PetscLimiterLimit_VanLeer;
768: return(0);
769: }
771: /*MC
772: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
774: Level: intermediate
776: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
777: M*/
781: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
782: {
783: PetscLimiter_VanLeer *l;
784: PetscErrorCode ierr;
788: PetscNewLog(lim, &l);
789: lim->data = l;
791: PetscLimiterInitialize_VanLeer(lim);
792: return(0);
793: }
797: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter fvm)
798: {
799: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) fvm->data;
800: PetscErrorCode ierr;
803: PetscFree(l);
804: return(0);
805: }
809: PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter fv, PetscViewer viewer)
810: {
811: PetscViewerFormat format;
812: PetscErrorCode ierr;
815: PetscViewerGetFormat(viewer, &format);
816: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
817: return(0);
818: }
822: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter fv, PetscViewer viewer)
823: {
824: PetscBool iascii;
830: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
831: if (iascii) {PetscLimiterView_VanAlbada_Ascii(fv, viewer);}
832: return(0);
833: }
837: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
838: {
840: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
841: return(0);
842: }
846: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter fvm)
847: {
849: fvm->ops->view = PetscLimiterView_VanAlbada;
850: fvm->ops->destroy = PetscLimiterDestroy_VanAlbada;
851: fvm->ops->limit = PetscLimiterLimit_VanAlbada;
852: return(0);
853: }
855: /*MC
856: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
858: Level: intermediate
860: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
861: M*/
865: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
866: {
867: PetscLimiter_VanAlbada *l;
868: PetscErrorCode ierr;
872: PetscNewLog(lim, &l);
873: lim->data = l;
875: PetscLimiterInitialize_VanAlbada(lim);
876: return(0);
877: }
881: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter fvm)
882: {
883: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) fvm->data;
884: PetscErrorCode ierr;
887: PetscFree(l);
888: return(0);
889: }
893: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter fv, PetscViewer viewer)
894: {
895: PetscViewerFormat format;
896: PetscErrorCode ierr;
899: PetscViewerGetFormat(viewer, &format);
900: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
901: return(0);
902: }
906: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter fv, PetscViewer viewer)
907: {
908: PetscBool iascii;
914: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
915: if (iascii) {PetscLimiterView_Superbee_Ascii(fv, viewer);}
916: return(0);
917: }
921: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
922: {
924: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
925: return(0);
926: }
930: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter fvm)
931: {
933: fvm->ops->view = PetscLimiterView_Superbee;
934: fvm->ops->destroy = PetscLimiterDestroy_Superbee;
935: fvm->ops->limit = PetscLimiterLimit_Superbee;
936: return(0);
937: }
939: /*MC
940: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
942: Level: intermediate
944: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
945: M*/
949: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
950: {
951: PetscLimiter_Superbee *l;
952: PetscErrorCode ierr;
956: PetscNewLog(lim, &l);
957: lim->data = l;
959: PetscLimiterInitialize_Superbee(lim);
960: return(0);
961: }
965: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter fvm)
966: {
967: PetscLimiter_MC *l = (PetscLimiter_MC *) fvm->data;
968: PetscErrorCode ierr;
971: PetscFree(l);
972: return(0);
973: }
977: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter fv, PetscViewer viewer)
978: {
979: PetscViewerFormat format;
980: PetscErrorCode ierr;
983: PetscViewerGetFormat(viewer, &format);
984: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
985: return(0);
986: }
990: PetscErrorCode PetscLimiterView_MC(PetscLimiter fv, PetscViewer viewer)
991: {
992: PetscBool iascii;
998: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
999: if (iascii) {PetscLimiterView_MC_Ascii(fv, viewer);}
1000: return(0);
1001: }
1005: /* aka Barth-Jespersen */
1006: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
1007: {
1009: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
1010: return(0);
1011: }
1015: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter fvm)
1016: {
1018: fvm->ops->view = PetscLimiterView_MC;
1019: fvm->ops->destroy = PetscLimiterDestroy_MC;
1020: fvm->ops->limit = PetscLimiterLimit_MC;
1021: return(0);
1022: }
1024: /*MC
1025: PETSCLIMITERMC = "mc" - A PetscLimiter object
1027: Level: intermediate
1029: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
1030: M*/
1034: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
1035: {
1036: PetscLimiter_MC *l;
1037: PetscErrorCode ierr;
1041: PetscNewLog(lim, &l);
1042: lim->data = l;
1044: PetscLimiterInitialize_MC(lim);
1045: return(0);
1046: }
1048: PetscClassId PETSCFV_CLASSID = 0;
1050: PetscFunctionList PetscFVList = NULL;
1051: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
1055: /*@C
1056: PetscFVRegister - Adds a new PetscFV implementation
1058: Not Collective
1060: Input Parameters:
1061: + name - The name of a new user-defined creation routine
1062: - create_func - The creation routine itself
1064: Notes:
1065: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
1067: Sample usage:
1068: .vb
1069: PetscFVRegister("my_fv", MyPetscFVCreate);
1070: .ve
1072: Then, your PetscFV type can be chosen with the procedural interface via
1073: .vb
1074: PetscFVCreate(MPI_Comm, PetscFV *);
1075: PetscFVSetType(PetscFV, "my_fv");
1076: .ve
1077: or at runtime via the option
1078: .vb
1079: -petscfv_type my_fv
1080: .ve
1082: Level: advanced
1084: .keywords: PetscFV, register
1085: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
1087: @*/
1088: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
1089: {
1093: PetscFunctionListAdd(&PetscFVList, sname, function);
1094: return(0);
1095: }
1099: /*@C
1100: PetscFVSetType - Builds a particular PetscFV
1102: Collective on PetscFV
1104: Input Parameters:
1105: + fvm - The PetscFV object
1106: - name - The kind of FVM space
1108: Options Database Key:
1109: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
1111: Level: intermediate
1113: .keywords: PetscFV, set, type
1114: .seealso: PetscFVGetType(), PetscFVCreate()
1115: @*/
1116: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
1117: {
1118: PetscErrorCode (*r)(PetscFV);
1119: PetscBool match;
1124: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1125: if (match) return(0);
1127: if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}
1128: PetscFunctionListFind(PetscFVList, name, &r);
1129: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1131: if (fvm->ops->destroy) {
1132: (*fvm->ops->destroy)(fvm);
1133: fvm->ops->destroy = NULL;
1134: }
1135: (*r)(fvm);
1136: PetscObjectChangeTypeName((PetscObject) fvm, name);
1137: return(0);
1138: }
1142: /*@C
1143: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1145: Not Collective
1147: Input Parameter:
1148: . fvm - The PetscFV
1150: Output Parameter:
1151: . name - The PetscFV type name
1153: Level: intermediate
1155: .keywords: PetscFV, get, type, name
1156: .seealso: PetscFVSetType(), PetscFVCreate()
1157: @*/
1158: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1159: {
1165: if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}
1166: *name = ((PetscObject) fvm)->type_name;
1167: return(0);
1168: }
1172: /*@C
1173: PetscFVView - Views a PetscFV
1175: Collective on PetscFV
1177: Input Parameter:
1178: + fvm - the PetscFV object to view
1179: - v - the viewer
1181: Level: developer
1183: .seealso: PetscFVDestroy()
1184: @*/
1185: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1186: {
1191: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1192: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1193: return(0);
1194: }
1198: /*
1199: PetscFVViewFromOptions - Processes command line options to determine if/how a PetscFV is to be viewed.
1201: Collective on PetscFV
1203: Input Parameters:
1204: + fvm - the PetscFV
1205: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
1206: - optionname - option to activate viewing
1208: Level: intermediate
1210: .keywords: PetscFV, view, options, database
1211: .seealso: VecViewFromOptions(), MatViewFromOptions()
1212: */
1213: PetscErrorCode PetscFVViewFromOptions(PetscFV fvm, const char prefix[], const char optionname[])
1214: {
1215: PetscViewer viewer;
1216: PetscViewerFormat format;
1217: PetscBool flg;
1218: PetscErrorCode ierr;
1221: if (prefix) {PetscOptionsGetViewer(PetscObjectComm((PetscObject) fvm), prefix, optionname, &viewer, &format, &flg);}
1222: else {PetscOptionsGetViewer(PetscObjectComm((PetscObject) fvm), ((PetscObject) fvm)->prefix, optionname, &viewer, &format, &flg);}
1223: if (flg) {
1224: PetscViewerPushFormat(viewer, format);
1225: PetscFVView(fvm, viewer);
1226: PetscViewerPopFormat(viewer);
1227: PetscViewerDestroy(&viewer);
1228: }
1229: return(0);
1230: }
1234: /*@
1235: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1237: Collective on PetscFV
1239: Input Parameter:
1240: . fvm - the PetscFV object to set options for
1242: Level: developer
1244: .seealso: PetscFVView()
1245: @*/
1246: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1247: {
1248: const char *defaultType;
1249: char name[256];
1250: PetscBool flg;
1255: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1256: else defaultType = ((PetscObject) fvm)->type_name;
1257: if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}
1259: PetscObjectOptionsBegin((PetscObject) fvm);
1260: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1261: if (flg) {
1262: PetscFVSetType(fvm, name);
1263: } else if (!((PetscObject) fvm)->type_name) {
1264: PetscFVSetType(fvm, defaultType);
1265: }
1266: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1267: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1268: PetscObjectProcessOptionsHandlers((PetscObject) fvm);
1269: PetscLimiterSetFromOptions(fvm->limiter);
1270: PetscOptionsEnd();
1271: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1272: return(0);
1273: }
1277: /*@
1278: PetscFVSetUp - Construct data structures for the PetscFV
1280: Collective on PetscFV
1282: Input Parameter:
1283: . fvm - the PetscFV object to setup
1285: Level: developer
1287: .seealso: PetscFVView(), PetscFVDestroy()
1288: @*/
1289: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1290: {
1295: PetscLimiterSetUp(fvm->limiter);
1296: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1297: return(0);
1298: }
1302: /*@
1303: PetscFVDestroy - Destroys a PetscFV object
1305: Collective on PetscFV
1307: Input Parameter:
1308: . fvm - the PetscFV object to destroy
1310: Level: developer
1312: .seealso: PetscFVView()
1313: @*/
1314: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1315: {
1319: if (!*fvm) return(0);
1322: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1323: ((PetscObject) (*fvm))->refct = 0;
1325: PetscLimiterDestroy(&(*fvm)->limiter);
1326: PetscFree((*fvm)->fluxWork);
1328: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1329: PetscHeaderDestroy(fvm);
1330: return(0);
1331: }
1335: /*@
1336: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1338: Collective on MPI_Comm
1340: Input Parameter:
1341: . comm - The communicator for the PetscFV object
1343: Output Parameter:
1344: . fvm - The PetscFV object
1346: Level: beginner
1348: .seealso: PetscFVSetType(), PETSCFVUPWIND
1349: @*/
1350: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1351: {
1352: PetscFV f;
1357: *fvm = NULL;
1358: PetscFVInitializePackage();
1360: PetscHeaderCreate(f, _p_PetscFV, struct _PetscFVOps, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1361: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1363: PetscLimiterCreate(comm, &f->limiter);
1364: f->numComponents = 1;
1365: f->dim = 0;
1366: f->computeGradients = PETSC_FALSE;
1367: f->fluxWork = NULL;
1369: *fvm = f;
1370: return(0);
1371: }
1375: /*@
1376: PetscFVSetLimiter - Set the limiter object
1378: Logically collective on PetscFV
1380: Input Parameters:
1381: + fvm - the PetscFV object to destroy
1382: - lim - The PetscLimiter
1384: Level: developer
1386: .seealso: PetscFVGetLimiter()
1387: @*/
1388: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1389: {
1395: PetscLimiterDestroy(&fvm->limiter);
1396: PetscObjectReference((PetscObject) lim);
1397: fvm->limiter = lim;
1398: return(0);
1399: }
1403: /*@
1404: PetscFVGetLimiter - Get the limiter object
1406: Not collective
1408: Input Parameter:
1409: . fvm - the PetscFV object to destroy
1411: Output Parameter:
1412: . lim - The PetscLimiter
1414: Level: developer
1416: .seealso: PetscFVSetLimiter()
1417: @*/
1418: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1419: {
1423: *lim = fvm->limiter;
1424: return(0);
1425: }
1429: /*@
1430: PetscFVSetNumComponents - Set the number of field components
1432: Logically collective on PetscFV
1434: Input Parameters:
1435: + fvm - the PetscFV object to destroy
1436: - comp - The number of components
1438: Level: developer
1440: .seealso: PetscFVGetNumComponents()
1441: @*/
1442: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1443: {
1448: fvm->numComponents = comp;
1449: PetscFree(fvm->fluxWork);
1450: PetscMalloc1(comp, &fvm->fluxWork);
1451: return(0);
1452: }
1456: /*@
1457: PetscFVGetNumComponents - Get the number of field components
1459: Not collective
1461: Input Parameter:
1462: . fvm - the PetscFV object to destroy
1464: Output Parameter:
1465: , comp - The number of components
1467: Level: developer
1469: .seealso: PetscFVSetNumComponents()
1470: @*/
1471: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1472: {
1476: *comp = fvm->numComponents;
1477: return(0);
1478: }
1482: /*@
1483: PetscFVSetSpatialDimension - Set the spatial dimension
1485: Logically collective on PetscFV
1487: Input Parameters:
1488: + fvm - the PetscFV object to destroy
1489: - dim - The spatial dimension
1491: Level: developer
1493: .seealso: PetscFVGetSpatialDimension()
1494: @*/
1495: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1496: {
1499: fvm->dim = dim;
1500: return(0);
1501: }
1505: /*@
1506: PetscFVGetSpatialDimension - Get the spatial dimension
1508: Logically collective on PetscFV
1510: Input Parameter:
1511: . fvm - the PetscFV object to destroy
1513: Output Parameter:
1514: . dim - The spatial dimension
1516: Level: developer
1518: .seealso: PetscFVSetSpatialDimension()
1519: @*/
1520: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1521: {
1525: *dim = fvm->dim;
1526: return(0);
1527: }
1531: /*@
1532: PetscFVSetComputeGradients - Toggle computation of cell gradients
1534: Logically collective on PetscFV
1536: Input Parameters:
1537: + fvm - the PetscFV object to destroy
1538: - computeGradients - Flag to compute cell gradients
1540: Level: developer
1542: .seealso: PetscFVGetComputeGradients()
1543: @*/
1544: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1545: {
1548: fvm->computeGradients = computeGradients;
1549: return(0);
1550: }
1554: /*@
1555: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1557: Not collective
1559: Input Parameter:
1560: . fvm - the PetscFV object to destroy
1562: Output Parameter:
1563: . computeGradients - Flag to compute cell gradients
1565: Level: developer
1567: .seealso: PetscFVSetComputeGradients()
1568: @*/
1569: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1570: {
1574: *computeGradients = fvm->computeGradients;
1575: return(0);
1576: }
1580: /*@C
1581: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1583: Input Parameters:
1584: + fvm - The PetscFV object
1585: . numFaces - The number of cell faces which are not constrained
1586: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1588: Level: developer
1590: .seealso: PetscFVCreate()
1591: @*/
1592: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1593: {
1598: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1599: return(0);
1600: }
1604: /*C
1605: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1607: Not collective
1609: Input Parameters:
1610: + fvm - The PetscFV object for the field being integrated
1611: . Nface - The number of faces in the chunk
1612: . Nf - The number of physical fields
1613: . fv - The PetscFV objects for each field
1614: . field - The field being integrated
1615: . fgeom - The face geometry for each face in the chunk
1616: . cgeom - The cell geometry for each pair of cells in the chunk
1617: . uL - The state from the cell on the left
1618: . uR - The state from the cell on the right
1619: . riemann - Riemann solver
1620: - ctx - User context passed to Riemann solve
1622: Output Parameter
1623: + fluxL - the left fluxes for each face
1624: - fluxR - the right fluxes for each face
1625: */
1626: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom, PetscScalar uL[], PetscScalar uR[],
1627: void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
1628: PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
1629: {
1634: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, Nfaces, Nf, fv, field, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, ctx);}
1635: return(0);
1636: }
1640: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1641: {
1642: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1646: PetscFree(b);
1647: return(0);
1648: }
1652: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1653: {
1654: PetscInt Nc;
1655: PetscViewerFormat format;
1656: PetscErrorCode ierr;
1659: PetscFVGetNumComponents(fv, &Nc);
1660: PetscViewerGetFormat(viewer, &format);
1661: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1662: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1663: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1664: } else {
1665: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1666: }
1667: return(0);
1668: }
1672: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1673: {
1674: PetscBool iascii;
1680: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1681: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1682: return(0);
1683: }
1687: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1688: {
1690: return(0);
1691: }
1695: /*
1696: fgeom[f]=>v0 is the centroid
1697: cgeom->vol[f*2+0] contains the left geom
1698: cgeom->vol[f*2+1] contains the right geom
1699: */
1700: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom,
1701: PetscScalar xL[], PetscScalar xR[],
1702: void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
1703: PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
1704: {
1705: PetscScalar *flux = fvm->fluxWork;
1706: PetscInt dim, pdim, f, d;
1710: PetscFVGetSpatialDimension(fvm, &dim);
1711: PetscFVGetNumComponents(fvm, &pdim);
1712: for (f = 0; f < Nfaces; ++f) {
1713: (*riemann)(&fgeom.v0[f*dim], &fgeom.n[f*dim], &xL[f*pdim], &xR[f*pdim], flux, ctx);
1714: for (d = 0; d < pdim; ++d) {
1715: fluxL[f*pdim+d] = flux[d] / cgeom.vol[f*2+0];
1716: fluxR[f*pdim+d] = flux[d] / cgeom.vol[f*2+1];
1717: }
1718: }
1719: return(0);
1720: }
1724: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1725: {
1727: fvm->ops->setfromoptions = NULL;
1728: fvm->ops->setup = PetscFVSetUp_Upwind;
1729: fvm->ops->view = PetscFVView_Upwind;
1730: fvm->ops->destroy = PetscFVDestroy_Upwind;
1731: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1732: return(0);
1733: }
1735: /*MC
1736: PETSCFVUPWIND = "upwind" - A PetscFV object
1738: Level: intermediate
1740: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1741: M*/
1745: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1746: {
1747: PetscFV_Upwind *b;
1748: PetscErrorCode ierr;
1752: PetscNewLog(fvm,&b);
1753: fvm->data = b;
1755: PetscFVInitialize_Upwind(fvm);
1756: return(0);
1757: }
1759: #include <petscblaslapack.h>
1763: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1764: {
1765: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1766: PetscErrorCode ierr;
1769: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1770: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1771: PetscFree(ls);
1772: return(0);
1773: }
1777: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1778: {
1779: PetscInt Nc;
1780: PetscViewerFormat format;
1781: PetscErrorCode ierr;
1784: PetscFVGetNumComponents(fv, &Nc);
1785: PetscViewerGetFormat(viewer, &format);
1786: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1787: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1788: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1789: } else {
1790: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1791: }
1792: return(0);
1793: }
1797: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1798: {
1799: PetscBool iascii;
1805: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1806: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1807: return(0);
1808: }
1812: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1813: {
1815: return(0);
1816: }
1820: /* Overwrites A. Can only handle full-rank problems with m>=n */
1821: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1822: {
1823: PetscBool debug = PETSC_FALSE;
1825: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1826: PetscScalar *R,*Q,*Aback,Alpha;
1829: if (debug) {
1830: PetscMalloc1(m*n,&Aback);
1831: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1832: }
1834: PetscBLASIntCast(m,&M);
1835: PetscBLASIntCast(n,&N);
1836: PetscBLASIntCast(mstride,&lda);
1837: PetscBLASIntCast(worksize,&ldwork);
1838: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1839: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1840: PetscFPTrapPop();
1841: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1842: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1844: /* Extract an explicit representation of Q */
1845: Q = Ainv;
1846: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1847: K = N; /* full rank */
1848: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
1849: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1851: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1852: Alpha = 1.0;
1853: ldb = lda;
1854: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1855: /* Ainv is Q, overwritten with inverse */
1857: if (debug) { /* Check that pseudo-inverse worked */
1858: PetscScalar Beta = 0.0;
1859: PetscBLASInt ldc;
1860: K = N;
1861: ldc = N;
1862: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1863: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1864: PetscFree(Aback);
1865: }
1866: return(0);
1867: }
1871: /* Overwrites A. Can handle degenerate problems and m<n. */
1872: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1873: {
1874: PetscBool debug = PETSC_FALSE;
1875: PetscScalar *Brhs, *Aback;
1876: PetscScalar *tmpwork;
1877: PetscReal rcond;
1878: PetscInt i, j, maxmn;
1879: PetscBLASInt M, N, nrhs, lda, ldb, irank, ldwork, info;
1883: if (debug) {
1884: PetscMalloc1(m*n,&Aback);
1885: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1886: }
1888: /* initialize to identity */
1889: tmpwork = Ainv;
1890: Brhs = work;
1891: maxmn = PetscMax(m,n);
1892: for (j=0; j<maxmn; j++) {
1893: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1894: }
1896: PetscBLASIntCast(m,&M);
1897: PetscBLASIntCast(n,&N);
1898: nrhs = M;
1899: PetscBLASIntCast(mstride,&lda);
1900: PetscBLASIntCast(maxmn,&ldb);
1901: PetscBLASIntCast(worksize,&ldwork);
1902: rcond = -1;
1903: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1904: #if defined(PETSC_USE_COMPLEX)
1905: if (tmpwork && rcond) rcond = 0.0; /* Get rid of compiler warning */
1906: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "I don't think this makes sense for complex numbers");
1907: #else
1908: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
1909: #endif
1910: PetscFPTrapPop();
1911: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
1912: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1913: 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");
1915: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
1916: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
1917: for (i=0; i<n; i++) {
1918: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
1919: }
1920: return(0);
1921: }
1923: #if 0
1926: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1927: {
1928: PetscReal grad[2] = {0, 0};
1929: const PetscInt *faces;
1930: PetscInt numFaces, f;
1931: PetscErrorCode ierr;
1934: DMPlexGetConeSize(dm, cell, &numFaces);
1935: DMPlexGetCone(dm, cell, &faces);
1936: for (f = 0; f < numFaces; ++f) {
1937: const PetscInt *fcells;
1938: const CellGeom *cg1;
1939: const FaceGeom *fg;
1941: DMPlexGetSupport(dm, faces[f], &fcells);
1942: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1943: for (i = 0; i < 2; ++i) {
1944: PetscScalar du;
1946: if (fcells[i] == c) continue;
1947: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1948: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1949: grad[0] += fg->grad[!i][0] * du;
1950: grad[1] += fg->grad[!i][1] * du;
1951: }
1952: }
1953: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1954: return(0);
1955: }
1956: #endif
1960: /*
1961: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1963: Input Parameters:
1964: + fvm - The PetscFV object
1965: . numFaces - The number of cell faces which are not constrained
1966: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1968: Level: developer
1970: .seealso: PetscFVCreate()
1971: */
1972: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1973: {
1974: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1975: const PetscBool useSVD = PETSC_TRUE;
1976: const PetscInt maxFaces = ls->maxFaces;
1977: PetscInt dim, f, d;
1978: PetscErrorCode ierr;
1981: if (numFaces > maxFaces) {
1982: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
1983: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
1984: }
1985: PetscFVGetSpatialDimension(fvm, &dim);
1986: for (f = 0; f < numFaces; ++f) {
1987: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
1988: }
1989: /* Overwrites B with garbage, returns Binv in row-major format */
1990: if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
1991: else {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
1992: for (f = 0; f < numFaces; ++f) {
1993: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
1994: }
1995: return(0);
1996: }
2000: /*
2001: fgeom[f]=>v0 is the centroid
2002: cgeom->vol[f*2+0] contains the left geom
2003: cgeom->vol[f*2+1] contains the right geom
2004: */
2005: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom,
2006: PetscScalar xL[], PetscScalar xR[],
2007: void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
2008: PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
2009: {
2010: PetscScalar *flux = fvm->fluxWork;
2011: PetscInt dim, pdim, f, d;
2015: PetscFVGetSpatialDimension(fvm, &dim);
2016: PetscFVGetNumComponents(fvm, &pdim);
2017: for (f = 0; f < Nfaces; ++f) {
2018: (*riemann)(&fgeom.v0[f*dim], &fgeom.n[f*dim], &xL[f*pdim], &xR[f*pdim], flux, ctx);
2019: for (d = 0; d < pdim; ++d) {
2020: fluxL[f*pdim+d] = flux[d] / cgeom.vol[f*2+0];
2021: fluxR[f*pdim+d] = flux[d] / cgeom.vol[f*2+1];
2022: }
2023: }
2024: return(0);
2025: }
2029: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2030: {
2031: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2032: PetscInt dim, m, n, nrhs, minwork;
2033: PetscErrorCode ierr;
2037: PetscFVGetSpatialDimension(fvm, &dim);
2038: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2039: ls->maxFaces = maxFaces;
2040: m = ls->maxFaces;
2041: n = dim;
2042: nrhs = ls->maxFaces;
2043: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2044: ls->workSize = 5*minwork; /* We can afford to be extra generous */
2045: PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2046: return(0);
2047: }
2051: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2052: {
2054: fvm->ops->setfromoptions = NULL;
2055: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2056: fvm->ops->view = PetscFVView_LeastSquares;
2057: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2058: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2059: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2060: return(0);
2061: }
2063: /*MC
2064: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2066: Level: intermediate
2068: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2069: M*/
2073: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2074: {
2075: PetscFV_LeastSquares *ls;
2076: PetscErrorCode ierr;
2080: PetscNewLog(fvm, &ls);
2081: fvm->data = ls;
2083: ls->maxFaces = -1;
2084: ls->workSize = -1;
2085: ls->B = NULL;
2086: ls->Binv = NULL;
2087: ls->tau = NULL;
2088: ls->work = NULL;
2090: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2091: PetscFVInitialize_LeastSquares(fvm);
2092: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2093: return(0);
2094: }
2098: /*@
2099: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2101: Not collective
2103: Input parameters:
2104: + fvm - The PetscFV object
2105: - maxFaces - The maximum number of cell faces
2107: Level: intermediate
2109: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2110: @*/
2111: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2112: {
2117: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2118: return(0);
2119: }