Actual source code: dtfv.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/petscfvimpl.h>
2: #include <petsc/private/dmpleximpl.h>
3: #include <petscds.h>
5: PetscClassId PETSCLIMITER_CLASSID = 0;
7: PetscFunctionList PetscLimiterList = NULL;
8: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
10: PetscBool Limitercite = PETSC_FALSE;
11: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
12: " title = {Analysis of slope limiters on irregular grids},\n"
13: " journal = {AIAA paper},\n"
14: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
15: " volume = {490},\n"
16: " year = {2005}\n}\n";
18: /*@C
19: PetscLimiterRegister - Adds a new PetscLimiter implementation
21: Not Collective
23: Input Parameters:
24: + name - The name of a new user-defined creation routine
25: - create_func - The creation routine itself
27: Notes:
28: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
30: Sample usage:
31: .vb
32: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
33: .ve
35: Then, your PetscLimiter type can be chosen with the procedural interface via
36: .vb
37: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
38: PetscLimiterSetType(PetscLimiter, "my_lim");
39: .ve
40: or at runtime via the option
41: .vb
42: -petsclimiter_type my_lim
43: .ve
45: Level: advanced
47: .keywords: PetscLimiter, register
48: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
50: @*/
51: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
52: {
56: PetscFunctionListAdd(&PetscLimiterList, sname, function);
57: return(0);
58: }
60: /*@C
61: PetscLimiterSetType - Builds a particular PetscLimiter
63: Collective on PetscLimiter
65: Input Parameters:
66: + lim - The PetscLimiter object
67: - name - The kind of limiter
69: Options Database Key:
70: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
72: Level: intermediate
74: .keywords: PetscLimiter, set, type
75: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
76: @*/
77: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
78: {
79: PetscErrorCode (*r)(PetscLimiter);
80: PetscBool match;
85: PetscObjectTypeCompare((PetscObject) lim, name, &match);
86: if (match) return(0);
88: PetscLimiterRegisterAll();
89: PetscFunctionListFind(PetscLimiterList, name, &r);
90: if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
92: if (lim->ops->destroy) {
93: (*lim->ops->destroy)(lim);
94: lim->ops->destroy = NULL;
95: }
96: (*r)(lim);
97: PetscObjectChangeTypeName((PetscObject) lim, name);
98: return(0);
99: }
101: /*@C
102: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
104: Not Collective
106: Input Parameter:
107: . lim - The PetscLimiter
109: Output Parameter:
110: . name - The PetscLimiter type name
112: Level: intermediate
114: .keywords: PetscLimiter, get, type, name
115: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
116: @*/
117: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
118: {
124: PetscLimiterRegisterAll();
125: *name = ((PetscObject) lim)->type_name;
126: return(0);
127: }
129: /*@C
130: PetscLimiterView - Views a PetscLimiter
132: Collective on PetscLimiter
134: Input Parameter:
135: + lim - the PetscLimiter object to view
136: - v - the viewer
138: Level: developer
140: .seealso: PetscLimiterDestroy()
141: @*/
142: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
143: {
148: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
149: if (lim->ops->view) {(*lim->ops->view)(lim, v);}
150: return(0);
151: }
153: /*@
154: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
156: Collective on PetscLimiter
158: Input Parameter:
159: . lim - the PetscLimiter object to set options for
161: Level: developer
163: .seealso: PetscLimiterView()
164: @*/
165: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
166: {
167: const char *defaultType;
168: char name[256];
169: PetscBool flg;
174: if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
175: else defaultType = ((PetscObject) lim)->type_name;
176: PetscLimiterRegisterAll();
178: PetscObjectOptionsBegin((PetscObject) lim);
179: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
180: if (flg) {
181: PetscLimiterSetType(lim, name);
182: } else if (!((PetscObject) lim)->type_name) {
183: PetscLimiterSetType(lim, defaultType);
184: }
185: if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
186: /* process any options handlers added with PetscObjectAddOptionsHandler() */
187: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
188: PetscOptionsEnd();
189: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
190: return(0);
191: }
193: /*@C
194: PetscLimiterSetUp - Construct data structures for the PetscLimiter
196: Collective on PetscLimiter
198: Input Parameter:
199: . lim - the PetscLimiter object to setup
201: Level: developer
203: .seealso: PetscLimiterView(), PetscLimiterDestroy()
204: @*/
205: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
206: {
211: if (lim->ops->setup) {(*lim->ops->setup)(lim);}
212: return(0);
213: }
215: /*@
216: PetscLimiterDestroy - Destroys a PetscLimiter object
218: Collective on PetscLimiter
220: Input Parameter:
221: . lim - the PetscLimiter object to destroy
223: Level: developer
225: .seealso: PetscLimiterView()
226: @*/
227: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
228: {
232: if (!*lim) return(0);
235: if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
236: ((PetscObject) (*lim))->refct = 0;
238: if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
239: PetscHeaderDestroy(lim);
240: return(0);
241: }
243: /*@
244: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
246: Collective on MPI_Comm
248: Input Parameter:
249: . comm - The communicator for the PetscLimiter object
251: Output Parameter:
252: . lim - The PetscLimiter object
254: Level: beginner
256: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
257: @*/
258: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
259: {
260: PetscLimiter l;
265: PetscCitationsRegister(LimiterCitation,&Limitercite);
266: *lim = NULL;
267: PetscFVInitializePackage();
269: PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
271: *lim = l;
272: return(0);
273: }
275: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
276: *
277: * The classical flux-limited formulation is psi(r) where
278: *
279: * r = (u[0] - u[-1]) / (u[1] - u[0])
280: *
281: * The second order TVD region is bounded by
282: *
283: * psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
284: *
285: * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
286: * phi(r)(r+1)/2 in which the reconstructed interface values are
287: *
288: * u(v) = u[0] + phi(r) (grad u)[0] v
289: *
290: * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
291: *
292: * 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))
293: *
294: * For a nicer symmetric formulation, rewrite in terms of
295: *
296: * f = (u[0] - u[-1]) / (u[1] - u[-1])
297: *
298: * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
299: *
300: * phi(r) = phi(1/r)
301: *
302: * becomes
303: *
304: * w(f) = w(1-f).
305: *
306: * The limiters below implement this final form w(f). The reference methods are
307: *
308: * w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
309: * */
310: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
311: {
317: (*lim->ops->limit)(lim, flim, phi);
318: return(0);
319: }
321: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
322: {
323: PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
324: PetscErrorCode ierr;
327: PetscFree(l);
328: return(0);
329: }
331: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
332: {
333: PetscViewerFormat format;
334: PetscErrorCode ierr;
337: PetscViewerGetFormat(viewer, &format);
338: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
339: return(0);
340: }
342: PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
343: {
344: PetscBool iascii;
350: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
351: if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
352: return(0);
353: }
355: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
356: {
358: *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
359: return(0);
360: }
362: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
363: {
365: lim->ops->view = PetscLimiterView_Sin;
366: lim->ops->destroy = PetscLimiterDestroy_Sin;
367: lim->ops->limit = PetscLimiterLimit_Sin;
368: return(0);
369: }
371: /*MC
372: PETSCLIMITERSIN = "sin" - A PetscLimiter object
374: Level: intermediate
376: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
377: M*/
379: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
380: {
381: PetscLimiter_Sin *l;
382: PetscErrorCode ierr;
386: PetscNewLog(lim, &l);
387: lim->data = l;
389: PetscLimiterInitialize_Sin(lim);
390: return(0);
391: }
393: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
394: {
395: PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
396: PetscErrorCode ierr;
399: PetscFree(l);
400: return(0);
401: }
403: PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
404: {
405: PetscViewerFormat format;
406: PetscErrorCode ierr;
409: PetscViewerGetFormat(viewer, &format);
410: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
411: return(0);
412: }
414: PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
415: {
416: PetscBool iascii;
422: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
423: if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
424: return(0);
425: }
427: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
428: {
430: *phi = 0.0;
431: return(0);
432: }
434: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
435: {
437: lim->ops->view = PetscLimiterView_Zero;
438: lim->ops->destroy = PetscLimiterDestroy_Zero;
439: lim->ops->limit = PetscLimiterLimit_Zero;
440: return(0);
441: }
443: /*MC
444: PETSCLIMITERZERO = "zero" - A PetscLimiter object
446: Level: intermediate
448: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
449: M*/
451: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
452: {
453: PetscLimiter_Zero *l;
454: PetscErrorCode ierr;
458: PetscNewLog(lim, &l);
459: lim->data = l;
461: PetscLimiterInitialize_Zero(lim);
462: return(0);
463: }
465: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
466: {
467: PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
468: PetscErrorCode ierr;
471: PetscFree(l);
472: return(0);
473: }
475: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
476: {
477: PetscViewerFormat format;
478: PetscErrorCode ierr;
481: PetscViewerGetFormat(viewer, &format);
482: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
483: return(0);
484: }
486: PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
487: {
488: PetscBool iascii;
494: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
495: if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
496: return(0);
497: }
499: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
500: {
502: *phi = 1.0;
503: return(0);
504: }
506: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
507: {
509: lim->ops->view = PetscLimiterView_None;
510: lim->ops->destroy = PetscLimiterDestroy_None;
511: lim->ops->limit = PetscLimiterLimit_None;
512: return(0);
513: }
515: /*MC
516: PETSCLIMITERNONE = "none" - A PetscLimiter object
518: Level: intermediate
520: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
521: M*/
523: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
524: {
525: PetscLimiter_None *l;
526: PetscErrorCode ierr;
530: PetscNewLog(lim, &l);
531: lim->data = l;
533: PetscLimiterInitialize_None(lim);
534: return(0);
535: }
537: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
538: {
539: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
540: PetscErrorCode ierr;
543: PetscFree(l);
544: return(0);
545: }
547: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
548: {
549: PetscViewerFormat format;
550: PetscErrorCode ierr;
553: PetscViewerGetFormat(viewer, &format);
554: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
555: return(0);
556: }
558: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
559: {
560: PetscBool iascii;
566: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
567: if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
568: return(0);
569: }
571: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
572: {
574: *phi = 2*PetscMax(0, PetscMin(f, 1-f));
575: return(0);
576: }
578: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
579: {
581: lim->ops->view = PetscLimiterView_Minmod;
582: lim->ops->destroy = PetscLimiterDestroy_Minmod;
583: lim->ops->limit = PetscLimiterLimit_Minmod;
584: return(0);
585: }
587: /*MC
588: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
590: Level: intermediate
592: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
593: M*/
595: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
596: {
597: PetscLimiter_Minmod *l;
598: PetscErrorCode ierr;
602: PetscNewLog(lim, &l);
603: lim->data = l;
605: PetscLimiterInitialize_Minmod(lim);
606: return(0);
607: }
609: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
610: {
611: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
612: PetscErrorCode ierr;
615: PetscFree(l);
616: return(0);
617: }
619: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
620: {
621: PetscViewerFormat format;
622: PetscErrorCode ierr;
625: PetscViewerGetFormat(viewer, &format);
626: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
627: return(0);
628: }
630: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
631: {
632: PetscBool iascii;
638: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
639: if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
640: return(0);
641: }
643: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
644: {
646: *phi = PetscMax(0, 4*f*(1-f));
647: return(0);
648: }
650: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
651: {
653: lim->ops->view = PetscLimiterView_VanLeer;
654: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
655: lim->ops->limit = PetscLimiterLimit_VanLeer;
656: return(0);
657: }
659: /*MC
660: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
662: Level: intermediate
664: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
665: M*/
667: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
668: {
669: PetscLimiter_VanLeer *l;
670: PetscErrorCode ierr;
674: PetscNewLog(lim, &l);
675: lim->data = l;
677: PetscLimiterInitialize_VanLeer(lim);
678: return(0);
679: }
681: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
682: {
683: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
684: PetscErrorCode ierr;
687: PetscFree(l);
688: return(0);
689: }
691: PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
692: {
693: PetscViewerFormat format;
694: PetscErrorCode ierr;
697: PetscViewerGetFormat(viewer, &format);
698: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
699: return(0);
700: }
702: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
703: {
704: PetscBool iascii;
710: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
711: if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
712: return(0);
713: }
715: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
716: {
718: *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
719: return(0);
720: }
722: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
723: {
725: lim->ops->view = PetscLimiterView_VanAlbada;
726: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
727: lim->ops->limit = PetscLimiterLimit_VanAlbada;
728: return(0);
729: }
731: /*MC
732: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
734: Level: intermediate
736: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
737: M*/
739: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
740: {
741: PetscLimiter_VanAlbada *l;
742: PetscErrorCode ierr;
746: PetscNewLog(lim, &l);
747: lim->data = l;
749: PetscLimiterInitialize_VanAlbada(lim);
750: return(0);
751: }
753: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
754: {
755: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
756: PetscErrorCode ierr;
759: PetscFree(l);
760: return(0);
761: }
763: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
764: {
765: PetscViewerFormat format;
766: PetscErrorCode ierr;
769: PetscViewerGetFormat(viewer, &format);
770: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
771: return(0);
772: }
774: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
775: {
776: PetscBool iascii;
782: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
783: if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
784: return(0);
785: }
787: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
788: {
790: *phi = 4*PetscMax(0, PetscMin(f, 1-f));
791: return(0);
792: }
794: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
795: {
797: lim->ops->view = PetscLimiterView_Superbee;
798: lim->ops->destroy = PetscLimiterDestroy_Superbee;
799: lim->ops->limit = PetscLimiterLimit_Superbee;
800: return(0);
801: }
803: /*MC
804: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
806: Level: intermediate
808: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
809: M*/
811: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
812: {
813: PetscLimiter_Superbee *l;
814: PetscErrorCode ierr;
818: PetscNewLog(lim, &l);
819: lim->data = l;
821: PetscLimiterInitialize_Superbee(lim);
822: return(0);
823: }
825: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
826: {
827: PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
828: PetscErrorCode ierr;
831: PetscFree(l);
832: return(0);
833: }
835: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
836: {
837: PetscViewerFormat format;
838: PetscErrorCode ierr;
841: PetscViewerGetFormat(viewer, &format);
842: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
843: return(0);
844: }
846: PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
847: {
848: PetscBool iascii;
854: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
855: if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
856: return(0);
857: }
859: /* aka Barth-Jespersen */
860: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
861: {
863: *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
864: return(0);
865: }
867: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
868: {
870: lim->ops->view = PetscLimiterView_MC;
871: lim->ops->destroy = PetscLimiterDestroy_MC;
872: lim->ops->limit = PetscLimiterLimit_MC;
873: return(0);
874: }
876: /*MC
877: PETSCLIMITERMC = "mc" - A PetscLimiter object
879: Level: intermediate
881: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
882: M*/
884: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
885: {
886: PetscLimiter_MC *l;
887: PetscErrorCode ierr;
891: PetscNewLog(lim, &l);
892: lim->data = l;
894: PetscLimiterInitialize_MC(lim);
895: return(0);
896: }
898: PetscClassId PETSCFV_CLASSID = 0;
900: PetscFunctionList PetscFVList = NULL;
901: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
903: /*@C
904: PetscFVRegister - Adds a new PetscFV implementation
906: Not Collective
908: Input Parameters:
909: + name - The name of a new user-defined creation routine
910: - create_func - The creation routine itself
912: Notes:
913: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
915: Sample usage:
916: .vb
917: PetscFVRegister("my_fv", MyPetscFVCreate);
918: .ve
920: Then, your PetscFV type can be chosen with the procedural interface via
921: .vb
922: PetscFVCreate(MPI_Comm, PetscFV *);
923: PetscFVSetType(PetscFV, "my_fv");
924: .ve
925: or at runtime via the option
926: .vb
927: -petscfv_type my_fv
928: .ve
930: Level: advanced
932: .keywords: PetscFV, register
933: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
935: @*/
936: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
937: {
941: PetscFunctionListAdd(&PetscFVList, sname, function);
942: return(0);
943: }
945: /*@C
946: PetscFVSetType - Builds a particular PetscFV
948: Collective on PetscFV
950: Input Parameters:
951: + fvm - The PetscFV object
952: - name - The kind of FVM space
954: Options Database Key:
955: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
957: Level: intermediate
959: .keywords: PetscFV, set, type
960: .seealso: PetscFVGetType(), PetscFVCreate()
961: @*/
962: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
963: {
964: PetscErrorCode (*r)(PetscFV);
965: PetscBool match;
970: PetscObjectTypeCompare((PetscObject) fvm, name, &match);
971: if (match) return(0);
973: PetscFVRegisterAll();
974: PetscFunctionListFind(PetscFVList, name, &r);
975: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
977: if (fvm->ops->destroy) {
978: (*fvm->ops->destroy)(fvm);
979: fvm->ops->destroy = NULL;
980: }
981: (*r)(fvm);
982: PetscObjectChangeTypeName((PetscObject) fvm, name);
983: return(0);
984: }
986: /*@C
987: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
989: Not Collective
991: Input Parameter:
992: . fvm - The PetscFV
994: Output Parameter:
995: . name - The PetscFV type name
997: Level: intermediate
999: .keywords: PetscFV, get, type, name
1000: .seealso: PetscFVSetType(), PetscFVCreate()
1001: @*/
1002: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1003: {
1009: PetscFVRegisterAll();
1010: *name = ((PetscObject) fvm)->type_name;
1011: return(0);
1012: }
1014: /*@C
1015: PetscFVView - Views a PetscFV
1017: Collective on PetscFV
1019: Input Parameter:
1020: + fvm - the PetscFV object to view
1021: - v - the viewer
1023: Level: developer
1025: .seealso: PetscFVDestroy()
1026: @*/
1027: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1028: {
1033: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1034: if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1035: return(0);
1036: }
1038: /*@
1039: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1041: Collective on PetscFV
1043: Input Parameter:
1044: . fvm - the PetscFV object to set options for
1046: Options Database Key:
1047: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1049: Level: developer
1051: .seealso: PetscFVView()
1052: @*/
1053: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1054: {
1055: const char *defaultType;
1056: char name[256];
1057: PetscBool flg;
1062: if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1063: else defaultType = ((PetscObject) fvm)->type_name;
1064: PetscFVRegisterAll();
1066: PetscObjectOptionsBegin((PetscObject) fvm);
1067: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1068: if (flg) {
1069: PetscFVSetType(fvm, name);
1070: } else if (!((PetscObject) fvm)->type_name) {
1071: PetscFVSetType(fvm, defaultType);
1073: }
1074: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1075: if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1076: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1077: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1078: PetscLimiterSetFromOptions(fvm->limiter);
1079: PetscOptionsEnd();
1080: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1081: return(0);
1082: }
1084: /*@
1085: PetscFVSetUp - Construct data structures for the PetscFV
1087: Collective on PetscFV
1089: Input Parameter:
1090: . fvm - the PetscFV object to setup
1092: Level: developer
1094: .seealso: PetscFVView(), PetscFVDestroy()
1095: @*/
1096: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1097: {
1102: PetscLimiterSetUp(fvm->limiter);
1103: if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1104: return(0);
1105: }
1107: /*@
1108: PetscFVDestroy - Destroys a PetscFV object
1110: Collective on PetscFV
1112: Input Parameter:
1113: . fvm - the PetscFV object to destroy
1115: Level: developer
1117: .seealso: PetscFVView()
1118: @*/
1119: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1120: {
1121: PetscInt i;
1125: if (!*fvm) return(0);
1128: if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1129: ((PetscObject) (*fvm))->refct = 0;
1131: for (i = 0; i < (*fvm)->numComponents; i++) {
1132: PetscFree((*fvm)->componentNames[i]);
1133: }
1134: PetscFree((*fvm)->componentNames);
1135: PetscLimiterDestroy(&(*fvm)->limiter);
1136: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1137: PetscFree((*fvm)->fluxWork);
1138: PetscQuadratureDestroy(&(*fvm)->quadrature);
1139: PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);
1141: if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1142: PetscHeaderDestroy(fvm);
1143: return(0);
1144: }
1146: /*@
1147: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1149: Collective on MPI_Comm
1151: Input Parameter:
1152: . comm - The communicator for the PetscFV object
1154: Output Parameter:
1155: . fvm - The PetscFV object
1157: Level: beginner
1159: .seealso: PetscFVSetType(), PETSCFVUPWIND
1160: @*/
1161: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1162: {
1163: PetscFV f;
1168: *fvm = NULL;
1169: PetscFVInitializePackage();
1171: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1172: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1174: PetscLimiterCreate(comm, &f->limiter);
1175: f->numComponents = 1;
1176: f->dim = 0;
1177: f->computeGradients = PETSC_FALSE;
1178: f->fluxWork = NULL;
1179: PetscCalloc1(f->numComponents,&f->componentNames);
1181: *fvm = f;
1182: return(0);
1183: }
1185: /*@
1186: PetscFVSetLimiter - Set the limiter object
1188: Logically collective on PetscFV
1190: Input Parameters:
1191: + fvm - the PetscFV object
1192: - lim - The PetscLimiter
1194: Level: developer
1196: .seealso: PetscFVGetLimiter()
1197: @*/
1198: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1199: {
1205: PetscLimiterDestroy(&fvm->limiter);
1206: PetscObjectReference((PetscObject) lim);
1207: fvm->limiter = lim;
1208: return(0);
1209: }
1211: /*@
1212: PetscFVGetLimiter - Get the limiter object
1214: Not collective
1216: Input Parameter:
1217: . fvm - the PetscFV object
1219: Output Parameter:
1220: . lim - The PetscLimiter
1222: Level: developer
1224: .seealso: PetscFVSetLimiter()
1225: @*/
1226: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1227: {
1231: *lim = fvm->limiter;
1232: return(0);
1233: }
1235: /*@
1236: PetscFVSetNumComponents - Set the number of field components
1238: Logically collective on PetscFV
1240: Input Parameters:
1241: + fvm - the PetscFV object
1242: - comp - The number of components
1244: Level: developer
1246: .seealso: PetscFVGetNumComponents()
1247: @*/
1248: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1249: {
1254: if (fvm->numComponents != comp) {
1255: PetscInt i;
1257: for (i = 0; i < fvm->numComponents; i++) {
1258: PetscFree(fvm->componentNames[i]);
1259: }
1260: PetscFree(fvm->componentNames);
1261: PetscCalloc1(comp,&fvm->componentNames);
1262: }
1263: fvm->numComponents = comp;
1264: PetscFree(fvm->fluxWork);
1265: PetscMalloc1(comp, &fvm->fluxWork);
1266: return(0);
1267: }
1269: /*@
1270: PetscFVGetNumComponents - Get the number of field components
1272: Not collective
1274: Input Parameter:
1275: . fvm - the PetscFV object
1277: Output Parameter:
1278: , comp - The number of components
1280: Level: developer
1282: .seealso: PetscFVSetNumComponents()
1283: @*/
1284: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1285: {
1289: *comp = fvm->numComponents;
1290: return(0);
1291: }
1293: /*@C
1294: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1296: Logically collective on PetscFV
1297: Input Parameters:
1298: + fvm - the PetscFV object
1299: . comp - the component number
1300: - name - the component name
1302: Level: developer
1304: .seealso: PetscFVGetComponentName()
1305: @*/
1306: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1307: {
1311: PetscFree(fvm->componentNames[comp]);
1312: PetscStrallocpy(name,&fvm->componentNames[comp]);
1313: return(0);
1314: }
1316: /*@C
1317: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1319: Logically collective on PetscFV
1320: Input Parameters:
1321: + fvm - the PetscFV object
1322: - comp - the component number
1324: Output Parameter:
1325: . name - the component name
1327: Level: developer
1329: .seealso: PetscFVSetComponentName()
1330: @*/
1331: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1332: {
1334: *name = fvm->componentNames[comp];
1335: return(0);
1336: }
1338: /*@
1339: PetscFVSetSpatialDimension - Set the spatial dimension
1341: Logically collective on PetscFV
1343: Input Parameters:
1344: + fvm - the PetscFV object
1345: - dim - The spatial dimension
1347: Level: developer
1349: .seealso: PetscFVGetSpatialDimension()
1350: @*/
1351: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1352: {
1355: fvm->dim = dim;
1356: return(0);
1357: }
1359: /*@
1360: PetscFVGetSpatialDimension - Get the spatial dimension
1362: Logically collective on PetscFV
1364: Input Parameter:
1365: . fvm - the PetscFV object
1367: Output Parameter:
1368: . dim - The spatial dimension
1370: Level: developer
1372: .seealso: PetscFVSetSpatialDimension()
1373: @*/
1374: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1375: {
1379: *dim = fvm->dim;
1380: return(0);
1381: }
1383: /*@
1384: PetscFVSetComputeGradients - Toggle computation of cell gradients
1386: Logically collective on PetscFV
1388: Input Parameters:
1389: + fvm - the PetscFV object
1390: - computeGradients - Flag to compute cell gradients
1392: Level: developer
1394: .seealso: PetscFVGetComputeGradients()
1395: @*/
1396: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1397: {
1400: fvm->computeGradients = computeGradients;
1401: return(0);
1402: }
1404: /*@
1405: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1407: Not collective
1409: Input Parameter:
1410: . fvm - the PetscFV object
1412: Output Parameter:
1413: . computeGradients - Flag to compute cell gradients
1415: Level: developer
1417: .seealso: PetscFVSetComputeGradients()
1418: @*/
1419: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1420: {
1424: *computeGradients = fvm->computeGradients;
1425: return(0);
1426: }
1428: /*@
1429: PetscFVSetQuadrature - Set the quadrature object
1431: Logically collective on PetscFV
1433: Input Parameters:
1434: + fvm - the PetscFV object
1435: - q - The PetscQuadrature
1437: Level: developer
1439: .seealso: PetscFVGetQuadrature()
1440: @*/
1441: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1442: {
1447: PetscQuadratureDestroy(&fvm->quadrature);
1448: PetscObjectReference((PetscObject) q);
1449: fvm->quadrature = q;
1450: return(0);
1451: }
1453: /*@
1454: PetscFVGetQuadrature - Get the quadrature object
1456: Not collective
1458: Input Parameter:
1459: . fvm - the PetscFV object
1461: Output Parameter:
1462: . lim - The PetscQuadrature
1464: Level: developer
1466: .seealso: PetscFVSetQuadrature()
1467: @*/
1468: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1469: {
1473: if (!fvm->quadrature) {
1474: /* Create default 1-point quadrature */
1475: PetscReal *points, *weights;
1478: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1479: PetscCalloc1(fvm->dim, &points);
1480: PetscMalloc1(1, &weights);
1481: weights[0] = 1.0;
1482: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1483: }
1484: *q = fvm->quadrature;
1485: return(0);
1486: }
1488: /*@
1489: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1491: Not collective
1493: Input Parameter:
1494: . fvm - The PetscFV object
1496: Output Parameter:
1497: . sp - The PetscDualSpace object
1499: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1501: Level: developer
1503: .seealso: PetscFVCreate()
1504: @*/
1505: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1506: {
1510: if (!fvm->dualSpace) {
1511: DM K;
1512: PetscInt dim, Nc, c;
1513: PetscErrorCode ierr;
1515: PetscFVGetSpatialDimension(fvm, &dim);
1516: PetscFVGetNumComponents(fvm, &Nc);
1517: PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1518: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1519: PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1520: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1521: PetscDualSpaceSetDM(fvm->dualSpace, K);
1522: DMDestroy(&K);
1523: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1524: /* Should we be using PetscFVGetQuadrature() here? */
1525: for (c = 0; c < Nc; ++c) {
1526: PetscQuadrature qc;
1527: PetscReal *points, *weights;
1528: PetscErrorCode ierr;
1530: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1531: PetscCalloc1(dim, &points);
1532: PetscCalloc1(Nc, &weights);
1533: weights[c] = 1.0;
1534: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1535: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1536: PetscQuadratureDestroy(&qc);
1537: }
1538: }
1539: *sp = fvm->dualSpace;
1540: return(0);
1541: }
1543: /*@
1544: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1546: Not collective
1548: Input Parameters:
1549: + fvm - The PetscFV object
1550: - sp - The PetscDualSpace object
1552: Level: intermediate
1554: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1556: .seealso: PetscFVCreate()
1557: @*/
1558: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1559: {
1565: PetscDualSpaceDestroy(&fvm->dualSpace);
1566: fvm->dualSpace = sp;
1567: PetscObjectReference((PetscObject) fvm->dualSpace);
1568: return(0);
1569: }
1571: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1572: {
1573: PetscInt npoints;
1574: const PetscReal *points;
1575: PetscErrorCode ierr;
1582: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1583: if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1584: if (B) *B = fvm->B;
1585: if (D) *D = fvm->D;
1586: if (H) *H = fvm->H;
1587: return(0);
1588: }
1590: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1591: {
1592: PetscInt pdim = 1; /* Dimension of approximation space P */
1593: PetscInt dim; /* Spatial dimension */
1594: PetscInt comp; /* Field components */
1595: PetscInt p, d, c, e;
1596: PetscErrorCode ierr;
1604: PetscFVGetSpatialDimension(fvm, &dim);
1605: PetscFVGetNumComponents(fvm, &comp);
1606: if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1607: if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1608: if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1609: 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;}
1610: 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;}
1611: 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;}
1612: return(0);
1613: }
1615: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1616: {
1621: if (B && *B) {PetscFree(*B);}
1622: if (D && *D) {PetscFree(*D);}
1623: if (H && *H) {PetscFree(*H);}
1624: return(0);
1625: }
1627: /*@C
1628: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1630: Input Parameters:
1631: + fvm - The PetscFV object
1632: . numFaces - The number of cell faces which are not constrained
1633: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
1635: Level: developer
1637: .seealso: PetscFVCreate()
1638: @*/
1639: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1640: {
1645: if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1646: return(0);
1647: }
1649: /*C
1650: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1652: Not collective
1654: Input Parameters:
1655: + fvm - The PetscFV object for the field being integrated
1656: . prob - The PetscDS specifing the discretizations and continuum functions
1657: . field - The field being integrated
1658: . Nf - The number of faces in the chunk
1659: . fgeom - The face geometry for each face in the chunk
1660: . neighborVol - The volume for each pair of cells in the chunk
1661: . uL - The state from the cell on the left
1662: - uR - The state from the cell on the right
1664: Output Parameter
1665: + fluxL - the left fluxes for each face
1666: - fluxR - the right fluxes for each face
1667: */
1668: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1669: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1670: {
1675: if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1676: return(0);
1677: }
1679: /*@
1680: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1681: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1682: sparsity). It is also used to create an interpolation between regularly refined meshes.
1684: Input Parameter:
1685: . fv - The initial PetscFV
1687: Output Parameter:
1688: . fvRef - The refined PetscFV
1690: Level: developer
1692: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1693: @*/
1694: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1695: {
1696: PetscDualSpace Q, Qref;
1697: DM K, Kref;
1698: PetscQuadrature q, qref;
1699: CellRefiner cellRefiner;
1700: PetscReal *v0;
1701: PetscReal *jac, *invjac;
1702: PetscInt numComp, numSubelements, s;
1703: PetscErrorCode ierr;
1706: PetscFVGetDualSpace(fv, &Q);
1707: PetscFVGetQuadrature(fv, &q);
1708: PetscDualSpaceGetDM(Q, &K);
1709: /* Create dual space */
1710: PetscDualSpaceDuplicate(Q, &Qref);
1711: DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1712: PetscDualSpaceSetDM(Qref, Kref);
1713: DMDestroy(&Kref);
1714: PetscDualSpaceSetUp(Qref);
1715: /* Create volume */
1716: PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1717: PetscFVSetDualSpace(*fvRef, Qref);
1718: PetscFVGetNumComponents(fv, &numComp);
1719: PetscFVSetNumComponents(*fvRef, numComp);
1720: PetscFVSetUp(*fvRef);
1721: /* Create quadrature */
1722: DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1723: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1724: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1725: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1726: for (s = 0; s < numSubelements; ++s) {
1727: PetscQuadrature qs;
1728: const PetscReal *points, *weights;
1729: PetscReal *p, *w;
1730: PetscInt dim, Nc, npoints, np;
1732: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1733: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1734: np = npoints/numSubelements;
1735: PetscMalloc1(np*dim,&p);
1736: PetscMalloc1(np*Nc,&w);
1737: PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1738: PetscMemcpy(w, &weights[s*np*Nc], np*Nc * sizeof(PetscReal));
1739: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1740: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1741: PetscQuadratureDestroy(&qs);
1742: }
1743: CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1744: PetscFVSetQuadrature(*fvRef, qref);
1745: PetscQuadratureDestroy(&qref);
1746: PetscDualSpaceDestroy(&Qref);
1747: return(0);
1748: }
1750: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1751: {
1752: PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1756: PetscFree(b);
1757: return(0);
1758: }
1760: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1761: {
1762: PetscInt Nc, c;
1763: PetscViewerFormat format;
1764: PetscErrorCode ierr;
1767: PetscFVGetNumComponents(fv, &Nc);
1768: PetscViewerGetFormat(viewer, &format);
1769: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1770: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1771: for (c = 0; c < Nc; c++) {
1772: if (fv->componentNames[c]) {
1773: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1774: }
1775: }
1776: return(0);
1777: }
1779: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1780: {
1781: PetscBool iascii;
1787: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1788: if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1789: return(0);
1790: }
1792: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1793: {
1795: return(0);
1796: }
1798: /*
1799: neighborVol[f*2+0] contains the left geom
1800: neighborVol[f*2+1] contains the right geom
1801: */
1802: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1803: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1804: {
1805: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1806: void *rctx;
1807: PetscScalar *flux = fvm->fluxWork;
1808: const PetscScalar *constants;
1809: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1810: PetscErrorCode ierr;
1813: PetscDSGetTotalComponents(prob, &Nc);
1814: PetscDSGetTotalDimension(prob, &totDim);
1815: PetscDSGetFieldOffset(prob, field, &off);
1816: PetscDSGetRiemannSolver(prob, field, &riemann);
1817: PetscDSGetContext(prob, field, &rctx);
1818: PetscDSGetConstants(prob, &numConstants, &constants);
1819: PetscFVGetSpatialDimension(fvm, &dim);
1820: PetscFVGetNumComponents(fvm, &pdim);
1821: for (f = 0; f < Nf; ++f) {
1822: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1823: for (d = 0; d < pdim; ++d) {
1824: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1825: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1826: }
1827: }
1828: return(0);
1829: }
1831: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1832: {
1834: fvm->ops->setfromoptions = NULL;
1835: fvm->ops->setup = PetscFVSetUp_Upwind;
1836: fvm->ops->view = PetscFVView_Upwind;
1837: fvm->ops->destroy = PetscFVDestroy_Upwind;
1838: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1839: return(0);
1840: }
1842: /*MC
1843: PETSCFVUPWIND = "upwind" - A PetscFV object
1845: Level: intermediate
1847: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1848: M*/
1850: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1851: {
1852: PetscFV_Upwind *b;
1853: PetscErrorCode ierr;
1857: PetscNewLog(fvm,&b);
1858: fvm->data = b;
1860: PetscFVInitialize_Upwind(fvm);
1861: return(0);
1862: }
1864: #include <petscblaslapack.h>
1866: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1867: {
1868: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1869: PetscErrorCode ierr;
1872: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1873: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1874: PetscFree(ls);
1875: return(0);
1876: }
1878: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1879: {
1880: PetscInt Nc, c;
1881: PetscViewerFormat format;
1882: PetscErrorCode ierr;
1885: PetscFVGetNumComponents(fv, &Nc);
1886: PetscViewerGetFormat(viewer, &format);
1887: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1888: PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);
1889: for (c = 0; c < Nc; c++) {
1890: if (fv->componentNames[c]) {
1891: PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);
1892: }
1893: }
1894: return(0);
1895: }
1897: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1898: {
1899: PetscBool iascii;
1905: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1906: if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1907: return(0);
1908: }
1910: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1911: {
1913: return(0);
1914: }
1916: /* Overwrites A. Can only handle full-rank problems with m>=n */
1917: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1918: {
1919: PetscBool debug = PETSC_FALSE;
1921: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1922: PetscScalar *R,*Q,*Aback,Alpha;
1925: if (debug) {
1926: PetscMalloc1(m*n,&Aback);
1927: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1928: }
1930: PetscBLASIntCast(m,&M);
1931: PetscBLASIntCast(n,&N);
1932: PetscBLASIntCast(mstride,&lda);
1933: PetscBLASIntCast(worksize,&ldwork);
1934: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1935: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1936: PetscFPTrapPop();
1937: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1938: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1940: /* Extract an explicit representation of Q */
1941: Q = Ainv;
1942: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1943: K = N; /* full rank */
1944: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
1945: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1947: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1948: Alpha = 1.0;
1949: ldb = lda;
1950: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1951: /* Ainv is Q, overwritten with inverse */
1953: if (debug) { /* Check that pseudo-inverse worked */
1954: PetscScalar Beta = 0.0;
1955: PetscBLASInt ldc;
1956: K = N;
1957: ldc = N;
1958: BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1959: PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1960: PetscFree(Aback);
1961: }
1962: return(0);
1963: }
1965: /* Overwrites A. Can handle degenerate problems and m<n. */
1966: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1967: {
1968: PetscBool debug = PETSC_FALSE;
1969: PetscScalar *Brhs, *Aback;
1970: PetscScalar *tmpwork;
1971: PetscReal rcond;
1972: #if defined (PETSC_USE_COMPLEX)
1973: PetscInt rworkSize;
1974: PetscReal *rwork;
1975: #endif
1976: PetscInt i, j, maxmn;
1977: PetscBLASInt M, N, lda, ldb, ldwork;
1978: PetscBLASInt nrhs, irank, info;
1982: if (debug) {
1983: PetscMalloc1(m*n,&Aback);
1984: PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1985: }
1987: /* initialize to identity */
1988: tmpwork = Ainv;
1989: Brhs = work;
1990: maxmn = PetscMax(m,n);
1991: for (j=0; j<maxmn; j++) {
1992: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1993: }
1995: PetscBLASIntCast(m,&M);
1996: PetscBLASIntCast(n,&N);
1997: PetscBLASIntCast(mstride,&lda);
1998: PetscBLASIntCast(maxmn,&ldb);
1999: PetscBLASIntCast(worksize,&ldwork);
2000: rcond = -1;
2001: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2002: nrhs = M;
2003: #if defined(PETSC_USE_COMPLEX)
2004: rworkSize = 5 * PetscMin(M,N);
2005: PetscMalloc1(rworkSize,&rwork);
2006: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2007: PetscFPTrapPop();
2008: PetscFree(rwork);
2009: #else
2010: nrhs = M;
2011: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2012: PetscFPTrapPop();
2013: #endif
2014: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2015: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2016: 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");
2018: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2019: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2020: for (i=0; i<n; i++) {
2021: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2022: }
2023: return(0);
2024: }
2026: #if 0
2027: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2028: {
2029: PetscReal grad[2] = {0, 0};
2030: const PetscInt *faces;
2031: PetscInt numFaces, f;
2032: PetscErrorCode ierr;
2035: DMPlexGetConeSize(dm, cell, &numFaces);
2036: DMPlexGetCone(dm, cell, &faces);
2037: for (f = 0; f < numFaces; ++f) {
2038: const PetscInt *fcells;
2039: const CellGeom *cg1;
2040: const FaceGeom *fg;
2042: DMPlexGetSupport(dm, faces[f], &fcells);
2043: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2044: for (i = 0; i < 2; ++i) {
2045: PetscScalar du;
2047: if (fcells[i] == c) continue;
2048: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2049: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2050: grad[0] += fg->grad[!i][0] * du;
2051: grad[1] += fg->grad[!i][1] * du;
2052: }
2053: }
2054: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2055: return(0);
2056: }
2057: #endif
2059: /*
2060: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2062: Input Parameters:
2063: + fvm - The PetscFV object
2064: . numFaces - The number of cell faces which are not constrained
2065: . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2067: Level: developer
2069: .seealso: PetscFVCreate()
2070: */
2071: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2072: {
2073: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2074: const PetscBool useSVD = PETSC_TRUE;
2075: const PetscInt maxFaces = ls->maxFaces;
2076: PetscInt dim, f, d;
2077: PetscErrorCode ierr;
2080: if (numFaces > maxFaces) {
2081: if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2082: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2083: }
2084: PetscFVGetSpatialDimension(fvm, &dim);
2085: for (f = 0; f < numFaces; ++f) {
2086: for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2087: }
2088: /* Overwrites B with garbage, returns Binv in row-major format */
2089: if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2090: else {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2091: for (f = 0; f < numFaces; ++f) {
2092: for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2093: }
2094: return(0);
2095: }
2097: /*
2098: neighborVol[f*2+0] contains the left geom
2099: neighborVol[f*2+1] contains the right geom
2100: */
2101: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2102: PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2103: {
2104: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2105: void *rctx;
2106: PetscScalar *flux = fvm->fluxWork;
2107: const PetscScalar *constants;
2108: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2109: PetscErrorCode ierr;
2112: PetscDSGetTotalComponents(prob, &Nc);
2113: PetscDSGetTotalDimension(prob, &totDim);
2114: PetscDSGetFieldOffset(prob, field, &off);
2115: PetscDSGetRiemannSolver(prob, field, &riemann);
2116: PetscDSGetContext(prob, field, &rctx);
2117: PetscDSGetConstants(prob, &numConstants, &constants);
2118: PetscFVGetSpatialDimension(fvm, &dim);
2119: PetscFVGetNumComponents(fvm, &pdim);
2120: for (f = 0; f < Nf; ++f) {
2121: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2122: for (d = 0; d < pdim; ++d) {
2123: fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2124: fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2125: }
2126: }
2127: return(0);
2128: }
2130: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2131: {
2132: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2133: PetscInt dim, m, n, nrhs, minwork;
2134: PetscErrorCode ierr;
2138: PetscFVGetSpatialDimension(fvm, &dim);
2139: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2140: ls->maxFaces = maxFaces;
2141: m = ls->maxFaces;
2142: n = dim;
2143: nrhs = ls->maxFaces;
2144: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2145: ls->workSize = 5*minwork; /* We can afford to be extra generous */
2146: PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2147: return(0);
2148: }
2150: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2151: {
2153: fvm->ops->setfromoptions = NULL;
2154: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2155: fvm->ops->view = PetscFVView_LeastSquares;
2156: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2157: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2158: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2159: return(0);
2160: }
2162: /*MC
2163: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2165: Level: intermediate
2167: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2168: M*/
2170: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2171: {
2172: PetscFV_LeastSquares *ls;
2173: PetscErrorCode ierr;
2177: PetscNewLog(fvm, &ls);
2178: fvm->data = ls;
2180: ls->maxFaces = -1;
2181: ls->workSize = -1;
2182: ls->B = NULL;
2183: ls->Binv = NULL;
2184: ls->tau = NULL;
2185: ls->work = NULL;
2187: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2188: PetscFVInitialize_LeastSquares(fvm);
2189: PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2190: return(0);
2191: }
2193: /*@
2194: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2196: Not collective
2198: Input parameters:
2199: + fvm - The PetscFV object
2200: - maxFaces - The maximum number of cell faces
2202: Level: intermediate
2204: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2205: @*/
2206: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2207: {
2212: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2213: return(0);
2214: }