Actual source code: fv.c
1: #include <petsc/private/petscfvimpl.h>
2: #include <petscdmplex.h>
3: #include <petscdmplextransform.h>
4: #include <petscds.h>
6: PetscClassId PETSCLIMITER_CLASSID = 0;
8: PetscFunctionList PetscLimiterList = NULL;
9: PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
11: PetscBool Limitercite = PETSC_FALSE;
12: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
13: " title = {Analysis of slope limiters on irregular grids},\n"
14: " journal = {AIAA paper},\n"
15: " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
16: " volume = {490},\n"
17: " year = {2005}\n}\n";
19: /*@C
20: PetscLimiterRegister - Adds a new PetscLimiter implementation
22: Not Collective
24: Input Parameters:
25: + name - The name of a new user-defined creation routine
26: - create_func - The creation routine itself
28: Notes:
29: PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
31: Sample usage:
32: .vb
33: PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
34: .ve
36: Then, your PetscLimiter type can be chosen with the procedural interface via
37: .vb
38: PetscLimiterCreate(MPI_Comm, PetscLimiter *);
39: PetscLimiterSetType(PetscLimiter, "my_lim");
40: .ve
41: or at runtime via the option
42: .vb
43: -petsclimiter_type my_lim
44: .ve
46: Level: advanced
48: .seealso: `PetscLimiterRegisterAll()`, `PetscLimiterRegisterDestroy()`
50: @*/
51: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
52: {
53: PetscFunctionListAdd(&PetscLimiterList, sname, function);
54: return 0;
55: }
57: /*@C
58: PetscLimiterSetType - Builds a particular PetscLimiter
60: Collective on lim
62: Input Parameters:
63: + lim - The PetscLimiter object
64: - name - The kind of limiter
66: Options Database Key:
67: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
69: Level: intermediate
71: .seealso: `PetscLimiterGetType()`, `PetscLimiterCreate()`
72: @*/
73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
74: {
75: PetscErrorCode (*r)(PetscLimiter);
76: PetscBool match;
79: PetscObjectTypeCompare((PetscObject)lim, name, &match);
80: if (match) return 0;
82: PetscLimiterRegisterAll();
83: PetscFunctionListFind(PetscLimiterList, name, &r);
86: if (lim->ops->destroy) {
87: PetscUseTypeMethod(lim, destroy);
88: lim->ops->destroy = NULL;
89: }
90: (*r)(lim);
91: PetscObjectChangeTypeName((PetscObject)lim, name);
92: return 0;
93: }
95: /*@C
96: PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
98: Not Collective
100: Input Parameter:
101: . lim - The PetscLimiter
103: Output Parameter:
104: . name - The PetscLimiter type name
106: Level: intermediate
108: .seealso: `PetscLimiterSetType()`, `PetscLimiterCreate()`
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
114: PetscLimiterRegisterAll();
115: *name = ((PetscObject)lim)->type_name;
116: return 0;
117: }
119: /*@C
120: PetscLimiterViewFromOptions - View from Options
122: Collective on PetscLimiter
124: Input Parameters:
125: + A - the PetscLimiter object to view
126: . obj - Optional object
127: - name - command line option
129: Level: intermediate
130: .seealso: `PetscLimiter`, `PetscLimiterView`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
131: @*/
132: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
133: {
135: PetscObjectViewFromOptions((PetscObject)A, obj, name);
136: return 0;
137: }
139: /*@C
140: PetscLimiterView - Views a PetscLimiter
142: Collective on lim
144: Input Parameters:
145: + lim - the PetscLimiter object to view
146: - v - the viewer
148: Level: beginner
150: .seealso: `PetscLimiterDestroy()`
151: @*/
152: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
153: {
155: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v);
156: PetscTryTypeMethod(lim, view, v);
157: return 0;
158: }
160: /*@
161: PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
163: Collective on lim
165: Input Parameter:
166: . lim - the PetscLimiter object to set options for
168: Level: intermediate
170: .seealso: `PetscLimiterView()`
171: @*/
172: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
173: {
174: const char *defaultType;
175: char name[256];
176: PetscBool flg;
179: if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
180: else defaultType = ((PetscObject)lim)->type_name;
181: PetscLimiterRegisterAll();
183: PetscObjectOptionsBegin((PetscObject)lim);
184: PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
185: if (flg) {
186: PetscLimiterSetType(lim, name);
187: } else if (!((PetscObject)lim)->type_name) {
188: PetscLimiterSetType(lim, defaultType);
189: }
190: PetscTryTypeMethod(lim, setfromoptions);
191: /* process any options handlers added with PetscObjectAddOptionsHandler() */
192: PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject);
193: PetscOptionsEnd();
194: PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
195: return 0;
196: }
198: /*@C
199: PetscLimiterSetUp - Construct data structures for the PetscLimiter
201: Collective on lim
203: Input Parameter:
204: . lim - the PetscLimiter object to setup
206: Level: intermediate
208: .seealso: `PetscLimiterView()`, `PetscLimiterDestroy()`
209: @*/
210: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
211: {
213: PetscTryTypeMethod(lim, setup);
214: return 0;
215: }
217: /*@
218: PetscLimiterDestroy - Destroys a PetscLimiter object
220: Collective on lim
222: Input Parameter:
223: . lim - the PetscLimiter object to destroy
225: Level: beginner
227: .seealso: `PetscLimiterView()`
228: @*/
229: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
230: {
231: if (!*lim) return 0;
234: if (--((PetscObject)(*lim))->refct > 0) {
235: *lim = NULL;
236: return 0;
237: }
238: ((PetscObject)(*lim))->refct = 0;
240: PetscTryTypeMethod((*lim), destroy);
241: PetscHeaderDestroy(lim);
242: return 0;
243: }
245: /*@
246: PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
248: Collective
250: Input Parameter:
251: . comm - The communicator for the PetscLimiter object
253: Output Parameter:
254: . lim - The PetscLimiter object
256: Level: beginner
258: .seealso: `PetscLimiterSetType()`, `PETSCLIMITERSIN`
259: @*/
260: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
261: {
262: 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: /*@
276: PetscLimiterLimit - Limit the flux
278: Input Parameters:
279: + lim - The PetscLimiter
280: - flim - The input field
282: Output Parameter:
283: . phi - The limited field
285: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
286: $ The classical flux-limited formulation is psi(r) where
287: $
288: $ r = (u[0] - u[-1]) / (u[1] - u[0])
289: $
290: $ The second order TVD region is bounded by
291: $
292: $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
293: $
294: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
295: $ phi(r)(r+1)/2 in which the reconstructed interface values are
296: $
297: $ u(v) = u[0] + phi(r) (grad u)[0] v
298: $
299: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
300: $
301: $ 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))
302: $
303: $ For a nicer symmetric formulation, rewrite in terms of
304: $
305: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
306: $
307: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
308: $
309: $ phi(r) = phi(1/r)
310: $
311: $ becomes
312: $
313: $ w(f) = w(1-f).
314: $
315: $ The limiters below implement this final form w(f). The reference methods are
316: $
317: $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
319: Level: beginner
321: .seealso: `PetscLimiterSetType()`, `PetscLimiterCreate()`
322: @*/
323: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
324: {
327: PetscUseTypeMethod(lim, limit, flim, phi);
328: return 0;
329: }
331: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
332: {
333: PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;
335: PetscFree(l);
336: return 0;
337: }
339: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
340: {
341: PetscViewerFormat format;
343: PetscViewerGetFormat(viewer, &format);
344: PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
345: return 0;
346: }
348: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
349: {
350: PetscBool iascii;
354: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
355: if (iascii) PetscLimiterView_Sin_Ascii(lim, viewer);
356: return 0;
357: }
359: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
360: {
361: *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
362: return 0;
363: }
365: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
366: {
367: lim->ops->view = PetscLimiterView_Sin;
368: lim->ops->destroy = PetscLimiterDestroy_Sin;
369: lim->ops->limit = PetscLimiterLimit_Sin;
370: return 0;
371: }
373: /*MC
374: PETSCLIMITERSIN = "sin" - A PetscLimiter object
376: Level: intermediate
378: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
379: M*/
381: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
382: {
383: PetscLimiter_Sin *l;
386: PetscNew(&l);
387: lim->data = l;
389: PetscLimiterInitialize_Sin(lim);
390: return 0;
391: }
393: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
394: {
395: PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;
397: PetscFree(l);
398: return 0;
399: }
401: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
402: {
403: PetscViewerFormat format;
405: PetscViewerGetFormat(viewer, &format);
406: PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
407: return 0;
408: }
410: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
411: {
412: PetscBool iascii;
416: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
417: if (iascii) PetscLimiterView_Zero_Ascii(lim, viewer);
418: return 0;
419: }
421: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
422: {
423: *phi = 0.0;
424: return 0;
425: }
427: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
428: {
429: lim->ops->view = PetscLimiterView_Zero;
430: lim->ops->destroy = PetscLimiterDestroy_Zero;
431: lim->ops->limit = PetscLimiterLimit_Zero;
432: return 0;
433: }
435: /*MC
436: PETSCLIMITERZERO = "zero" - A PetscLimiter object
438: Level: intermediate
440: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
441: M*/
443: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
444: {
445: PetscLimiter_Zero *l;
448: PetscNew(&l);
449: lim->data = l;
451: PetscLimiterInitialize_Zero(lim);
452: return 0;
453: }
455: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
456: {
457: PetscLimiter_None *l = (PetscLimiter_None *)lim->data;
459: PetscFree(l);
460: return 0;
461: }
463: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
464: {
465: PetscViewerFormat format;
467: PetscViewerGetFormat(viewer, &format);
468: PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
469: return 0;
470: }
472: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
473: {
474: PetscBool iascii;
478: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
479: if (iascii) PetscLimiterView_None_Ascii(lim, viewer);
480: return 0;
481: }
483: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
484: {
485: *phi = 1.0;
486: return 0;
487: }
489: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
490: {
491: lim->ops->view = PetscLimiterView_None;
492: lim->ops->destroy = PetscLimiterDestroy_None;
493: lim->ops->limit = PetscLimiterLimit_None;
494: return 0;
495: }
497: /*MC
498: PETSCLIMITERNONE = "none" - A PetscLimiter object
500: Level: intermediate
502: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
503: M*/
505: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
506: {
507: PetscLimiter_None *l;
510: PetscNew(&l);
511: lim->data = l;
513: PetscLimiterInitialize_None(lim);
514: return 0;
515: }
517: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
518: {
519: PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;
521: PetscFree(l);
522: return 0;
523: }
525: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
526: {
527: PetscViewerFormat format;
529: PetscViewerGetFormat(viewer, &format);
530: PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
531: return 0;
532: }
534: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
535: {
536: PetscBool iascii;
540: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
541: if (iascii) PetscLimiterView_Minmod_Ascii(lim, viewer);
542: return 0;
543: }
545: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
546: {
547: *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
548: return 0;
549: }
551: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
552: {
553: lim->ops->view = PetscLimiterView_Minmod;
554: lim->ops->destroy = PetscLimiterDestroy_Minmod;
555: lim->ops->limit = PetscLimiterLimit_Minmod;
556: return 0;
557: }
559: /*MC
560: PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
562: Level: intermediate
564: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
565: M*/
567: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
568: {
569: PetscLimiter_Minmod *l;
572: PetscNew(&l);
573: lim->data = l;
575: PetscLimiterInitialize_Minmod(lim);
576: return 0;
577: }
579: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
580: {
581: PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;
583: PetscFree(l);
584: return 0;
585: }
587: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
588: {
589: PetscViewerFormat format;
591: PetscViewerGetFormat(viewer, &format);
592: PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
593: return 0;
594: }
596: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
597: {
598: PetscBool iascii;
602: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
603: if (iascii) PetscLimiterView_VanLeer_Ascii(lim, viewer);
604: return 0;
605: }
607: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
608: {
609: *phi = PetscMax(0, 4 * f * (1 - f));
610: return 0;
611: }
613: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
614: {
615: lim->ops->view = PetscLimiterView_VanLeer;
616: lim->ops->destroy = PetscLimiterDestroy_VanLeer;
617: lim->ops->limit = PetscLimiterLimit_VanLeer;
618: return 0;
619: }
621: /*MC
622: PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
624: Level: intermediate
626: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
627: M*/
629: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
630: {
631: PetscLimiter_VanLeer *l;
634: PetscNew(&l);
635: lim->data = l;
637: PetscLimiterInitialize_VanLeer(lim);
638: return 0;
639: }
641: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
642: {
643: PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;
645: PetscFree(l);
646: return 0;
647: }
649: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
650: {
651: PetscViewerFormat format;
653: PetscViewerGetFormat(viewer, &format);
654: PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
655: return 0;
656: }
658: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
659: {
660: PetscBool iascii;
664: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
665: if (iascii) PetscLimiterView_VanAlbada_Ascii(lim, viewer);
666: return 0;
667: }
669: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
670: {
671: *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
672: return 0;
673: }
675: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
676: {
677: lim->ops->view = PetscLimiterView_VanAlbada;
678: lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
679: lim->ops->limit = PetscLimiterLimit_VanAlbada;
680: return 0;
681: }
683: /*MC
684: PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
686: Level: intermediate
688: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
689: M*/
691: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
692: {
693: PetscLimiter_VanAlbada *l;
696: PetscNew(&l);
697: lim->data = l;
699: PetscLimiterInitialize_VanAlbada(lim);
700: return 0;
701: }
703: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
704: {
705: PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;
707: PetscFree(l);
708: return 0;
709: }
711: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
712: {
713: PetscViewerFormat format;
715: PetscViewerGetFormat(viewer, &format);
716: PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
717: return 0;
718: }
720: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
721: {
722: PetscBool iascii;
726: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
727: if (iascii) PetscLimiterView_Superbee_Ascii(lim, viewer);
728: return 0;
729: }
731: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
732: {
733: *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
734: return 0;
735: }
737: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
738: {
739: lim->ops->view = PetscLimiterView_Superbee;
740: lim->ops->destroy = PetscLimiterDestroy_Superbee;
741: lim->ops->limit = PetscLimiterLimit_Superbee;
742: return 0;
743: }
745: /*MC
746: PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
748: Level: intermediate
750: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
751: M*/
753: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
754: {
755: PetscLimiter_Superbee *l;
758: PetscNew(&l);
759: lim->data = l;
761: PetscLimiterInitialize_Superbee(lim);
762: return 0;
763: }
765: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
766: {
767: PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;
769: PetscFree(l);
770: return 0;
771: }
773: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
774: {
775: PetscViewerFormat format;
777: PetscViewerGetFormat(viewer, &format);
778: PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
779: return 0;
780: }
782: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
783: {
784: PetscBool iascii;
788: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
789: if (iascii) PetscLimiterView_MC_Ascii(lim, viewer);
790: return 0;
791: }
793: /* aka Barth-Jespersen */
794: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
795: {
796: *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
797: return 0;
798: }
800: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
801: {
802: lim->ops->view = PetscLimiterView_MC;
803: lim->ops->destroy = PetscLimiterDestroy_MC;
804: lim->ops->limit = PetscLimiterLimit_MC;
805: return 0;
806: }
808: /*MC
809: PETSCLIMITERMC = "mc" - A PetscLimiter object
811: Level: intermediate
813: .seealso: `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
814: M*/
816: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
817: {
818: PetscLimiter_MC *l;
821: PetscNew(&l);
822: lim->data = l;
824: PetscLimiterInitialize_MC(lim);
825: return 0;
826: }
828: PetscClassId PETSCFV_CLASSID = 0;
830: PetscFunctionList PetscFVList = NULL;
831: PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
833: /*@C
834: PetscFVRegister - Adds a new PetscFV implementation
836: Not Collective
838: Input Parameters:
839: + name - The name of a new user-defined creation routine
840: - create_func - The creation routine itself
842: Notes:
843: PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
845: Sample usage:
846: .vb
847: PetscFVRegister("my_fv", MyPetscFVCreate);
848: .ve
850: Then, your PetscFV type can be chosen with the procedural interface via
851: .vb
852: PetscFVCreate(MPI_Comm, PetscFV *);
853: PetscFVSetType(PetscFV, "my_fv");
854: .ve
855: or at runtime via the option
856: .vb
857: -petscfv_type my_fv
858: .ve
860: Level: advanced
862: .seealso: `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
864: @*/
865: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
866: {
867: PetscFunctionListAdd(&PetscFVList, sname, function);
868: return 0;
869: }
871: /*@C
872: PetscFVSetType - Builds a particular PetscFV
874: Collective on fvm
876: Input Parameters:
877: + fvm - The PetscFV object
878: - name - The kind of FVM space
880: Options Database Key:
881: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
883: Level: intermediate
885: .seealso: `PetscFVGetType()`, `PetscFVCreate()`
886: @*/
887: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
888: {
889: PetscErrorCode (*r)(PetscFV);
890: PetscBool match;
893: PetscObjectTypeCompare((PetscObject)fvm, name, &match);
894: if (match) return 0;
896: PetscFVRegisterAll();
897: PetscFunctionListFind(PetscFVList, name, &r);
900: PetscTryTypeMethod(fvm, destroy);
901: fvm->ops->destroy = NULL;
903: (*r)(fvm);
904: PetscObjectChangeTypeName((PetscObject)fvm, name);
905: return 0;
906: }
908: /*@C
909: PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
911: Not Collective
913: Input Parameter:
914: . fvm - The PetscFV
916: Output Parameter:
917: . name - The PetscFV type name
919: Level: intermediate
921: .seealso: `PetscFVSetType()`, `PetscFVCreate()`
922: @*/
923: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
924: {
927: PetscFVRegisterAll();
928: *name = ((PetscObject)fvm)->type_name;
929: return 0;
930: }
932: /*@C
933: PetscFVViewFromOptions - View from Options
935: Collective on PetscFV
937: Input Parameters:
938: + A - the PetscFV object
939: . obj - Optional object
940: - name - command line option
942: Level: intermediate
943: .seealso: `PetscFV`, `PetscFVView`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
944: @*/
945: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
946: {
948: PetscObjectViewFromOptions((PetscObject)A, obj, name);
949: return 0;
950: }
952: /*@C
953: PetscFVView - Views a PetscFV
955: Collective on fvm
957: Input Parameters:
958: + fvm - the PetscFV object to view
959: - v - the viewer
961: Level: beginner
963: .seealso: `PetscFVDestroy()`
964: @*/
965: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
966: {
968: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v);
969: PetscTryTypeMethod(fvm, view, v);
970: return 0;
971: }
973: /*@
974: PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
976: Collective on fvm
978: Input Parameter:
979: . fvm - the PetscFV object to set options for
981: Options Database Key:
982: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
984: Level: intermediate
986: .seealso: `PetscFVView()`
987: @*/
988: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
989: {
990: const char *defaultType;
991: char name[256];
992: PetscBool flg;
995: if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
996: else defaultType = ((PetscObject)fvm)->type_name;
997: PetscFVRegisterAll();
999: PetscObjectOptionsBegin((PetscObject)fvm);
1000: PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1001: if (flg) {
1002: PetscFVSetType(fvm, name);
1003: } else if (!((PetscObject)fvm)->type_name) {
1004: PetscFVSetType(fvm, defaultType);
1005: }
1006: PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1007: PetscTryTypeMethod(fvm, setfromoptions);
1008: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1009: PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject);
1010: PetscLimiterSetFromOptions(fvm->limiter);
1011: PetscOptionsEnd();
1012: PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1013: return 0;
1014: }
1016: /*@
1017: PetscFVSetUp - Construct data structures for the PetscFV
1019: Collective on fvm
1021: Input Parameter:
1022: . fvm - the PetscFV object to setup
1024: Level: intermediate
1026: .seealso: `PetscFVView()`, `PetscFVDestroy()`
1027: @*/
1028: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1029: {
1031: PetscLimiterSetUp(fvm->limiter);
1032: PetscTryTypeMethod(fvm, setup);
1033: return 0;
1034: }
1036: /*@
1037: PetscFVDestroy - Destroys a PetscFV object
1039: Collective on fvm
1041: Input Parameter:
1042: . fvm - the PetscFV object to destroy
1044: Level: beginner
1046: .seealso: `PetscFVView()`
1047: @*/
1048: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1049: {
1050: PetscInt i;
1052: if (!*fvm) return 0;
1055: if (--((PetscObject)(*fvm))->refct > 0) {
1056: *fvm = NULL;
1057: return 0;
1058: }
1059: ((PetscObject)(*fvm))->refct = 0;
1061: for (i = 0; i < (*fvm)->numComponents; i++) PetscFree((*fvm)->componentNames[i]);
1062: PetscFree((*fvm)->componentNames);
1063: PetscLimiterDestroy(&(*fvm)->limiter);
1064: PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1065: PetscFree((*fvm)->fluxWork);
1066: PetscQuadratureDestroy(&(*fvm)->quadrature);
1067: PetscTabulationDestroy(&(*fvm)->T);
1069: PetscTryTypeMethod((*fvm), destroy);
1070: PetscHeaderDestroy(fvm);
1071: return 0;
1072: }
1074: /*@
1075: PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1077: Collective
1079: Input Parameter:
1080: . comm - The communicator for the PetscFV object
1082: Output Parameter:
1083: . fvm - The PetscFV object
1085: Level: beginner
1087: .seealso: `PetscFVSetType()`, `PETSCFVUPWIND`
1088: @*/
1089: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1090: {
1091: PetscFV f;
1094: *fvm = NULL;
1095: PetscFVInitializePackage();
1097: PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);
1098: PetscMemzero(f->ops, sizeof(struct _PetscFVOps));
1100: PetscLimiterCreate(comm, &f->limiter);
1101: f->numComponents = 1;
1102: f->dim = 0;
1103: f->computeGradients = PETSC_FALSE;
1104: f->fluxWork = NULL;
1105: PetscCalloc1(f->numComponents, &f->componentNames);
1107: *fvm = f;
1108: return 0;
1109: }
1111: /*@
1112: PetscFVSetLimiter - Set the limiter object
1114: Logically collective on fvm
1116: Input Parameters:
1117: + fvm - the PetscFV object
1118: - lim - The PetscLimiter
1120: Level: intermediate
1122: .seealso: `PetscFVGetLimiter()`
1123: @*/
1124: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1125: {
1128: PetscLimiterDestroy(&fvm->limiter);
1129: PetscObjectReference((PetscObject)lim);
1130: fvm->limiter = lim;
1131: return 0;
1132: }
1134: /*@
1135: PetscFVGetLimiter - Get the limiter object
1137: Not collective
1139: Input Parameter:
1140: . fvm - the PetscFV object
1142: Output Parameter:
1143: . lim - The PetscLimiter
1145: Level: intermediate
1147: .seealso: `PetscFVSetLimiter()`
1148: @*/
1149: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1150: {
1153: *lim = fvm->limiter;
1154: return 0;
1155: }
1157: /*@
1158: PetscFVSetNumComponents - Set the number of field components
1160: Logically collective on fvm
1162: Input Parameters:
1163: + fvm - the PetscFV object
1164: - comp - The number of components
1166: Level: intermediate
1168: .seealso: `PetscFVGetNumComponents()`
1169: @*/
1170: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1171: {
1173: if (fvm->numComponents != comp) {
1174: PetscInt i;
1176: for (i = 0; i < fvm->numComponents; i++) PetscFree(fvm->componentNames[i]);
1177: PetscFree(fvm->componentNames);
1178: PetscCalloc1(comp, &fvm->componentNames);
1179: }
1180: fvm->numComponents = comp;
1181: PetscFree(fvm->fluxWork);
1182: PetscMalloc1(comp, &fvm->fluxWork);
1183: return 0;
1184: }
1186: /*@
1187: PetscFVGetNumComponents - Get the number of field components
1189: Not collective
1191: Input Parameter:
1192: . fvm - the PetscFV object
1194: Output Parameter:
1195: , comp - The number of components
1197: Level: intermediate
1199: .seealso: `PetscFVSetNumComponents()`
1200: @*/
1201: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1202: {
1205: *comp = fvm->numComponents;
1206: return 0;
1207: }
1209: /*@C
1210: PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1212: Logically collective on fvm
1213: Input Parameters:
1214: + fvm - the PetscFV object
1215: . comp - the component number
1216: - name - the component name
1218: Level: intermediate
1220: .seealso: `PetscFVGetComponentName()`
1221: @*/
1222: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1223: {
1224: PetscFree(fvm->componentNames[comp]);
1225: PetscStrallocpy(name, &fvm->componentNames[comp]);
1226: return 0;
1227: }
1229: /*@C
1230: PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1232: Logically collective on fvm
1233: Input Parameters:
1234: + fvm - the PetscFV object
1235: - comp - the component number
1237: Output Parameter:
1238: . name - the component name
1240: Level: intermediate
1242: .seealso: `PetscFVSetComponentName()`
1243: @*/
1244: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1245: {
1246: *name = fvm->componentNames[comp];
1247: return 0;
1248: }
1250: /*@
1251: PetscFVSetSpatialDimension - Set the spatial dimension
1253: Logically collective on fvm
1255: Input Parameters:
1256: + fvm - the PetscFV object
1257: - dim - The spatial dimension
1259: Level: intermediate
1261: .seealso: `PetscFVGetSpatialDimension()`
1262: @*/
1263: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1264: {
1266: fvm->dim = dim;
1267: return 0;
1268: }
1270: /*@
1271: PetscFVGetSpatialDimension - Get the spatial dimension
1273: Logically collective on fvm
1275: Input Parameter:
1276: . fvm - the PetscFV object
1278: Output Parameter:
1279: . dim - The spatial dimension
1281: Level: intermediate
1283: .seealso: `PetscFVSetSpatialDimension()`
1284: @*/
1285: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1286: {
1289: *dim = fvm->dim;
1290: return 0;
1291: }
1293: /*@
1294: PetscFVSetComputeGradients - Toggle computation of cell gradients
1296: Logically collective on fvm
1298: Input Parameters:
1299: + fvm - the PetscFV object
1300: - computeGradients - Flag to compute cell gradients
1302: Level: intermediate
1304: .seealso: `PetscFVGetComputeGradients()`
1305: @*/
1306: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1307: {
1309: fvm->computeGradients = computeGradients;
1310: return 0;
1311: }
1313: /*@
1314: PetscFVGetComputeGradients - Return flag for computation of cell gradients
1316: Not collective
1318: Input Parameter:
1319: . fvm - the PetscFV object
1321: Output Parameter:
1322: . computeGradients - Flag to compute cell gradients
1324: Level: intermediate
1326: .seealso: `PetscFVSetComputeGradients()`
1327: @*/
1328: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1329: {
1332: *computeGradients = fvm->computeGradients;
1333: return 0;
1334: }
1336: /*@
1337: PetscFVSetQuadrature - Set the quadrature object
1339: Logically collective on fvm
1341: Input Parameters:
1342: + fvm - the PetscFV object
1343: - q - The PetscQuadrature
1345: Level: intermediate
1347: .seealso: `PetscFVGetQuadrature()`
1348: @*/
1349: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1350: {
1352: PetscQuadratureDestroy(&fvm->quadrature);
1353: PetscObjectReference((PetscObject)q);
1354: fvm->quadrature = q;
1355: return 0;
1356: }
1358: /*@
1359: PetscFVGetQuadrature - Get the quadrature object
1361: Not collective
1363: Input Parameter:
1364: . fvm - the PetscFV object
1366: Output Parameter:
1367: . lim - The PetscQuadrature
1369: Level: intermediate
1371: .seealso: `PetscFVSetQuadrature()`
1372: @*/
1373: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1374: {
1377: if (!fvm->quadrature) {
1378: /* Create default 1-point quadrature */
1379: PetscReal *points, *weights;
1381: PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1382: PetscCalloc1(fvm->dim, &points);
1383: PetscMalloc1(1, &weights);
1384: weights[0] = 1.0;
1385: PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1386: }
1387: *q = fvm->quadrature;
1388: return 0;
1389: }
1391: /*@
1392: PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1394: Not collective
1396: Input Parameter:
1397: . fvm - The PetscFV object
1399: Output Parameter:
1400: . sp - The PetscDualSpace object
1402: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1404: Level: intermediate
1406: .seealso: `PetscFVCreate()`
1407: @*/
1408: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1409: {
1412: if (!fvm->dualSpace) {
1413: DM K;
1414: PetscInt dim, Nc, c;
1416: PetscFVGetSpatialDimension(fvm, &dim);
1417: PetscFVGetNumComponents(fvm, &Nc);
1418: PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace);
1419: PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1420: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K);
1421: PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1422: PetscDualSpaceSetDM(fvm->dualSpace, K);
1423: DMDestroy(&K);
1424: PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1425: /* Should we be using PetscFVGetQuadrature() here? */
1426: for (c = 0; c < Nc; ++c) {
1427: PetscQuadrature qc;
1428: PetscReal *points, *weights;
1430: PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1431: PetscCalloc1(dim, &points);
1432: PetscCalloc1(Nc, &weights);
1433: weights[c] = 1.0;
1434: PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1435: PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1436: PetscQuadratureDestroy(&qc);
1437: }
1438: PetscDualSpaceSetUp(fvm->dualSpace);
1439: }
1440: *sp = fvm->dualSpace;
1441: return 0;
1442: }
1444: /*@
1445: PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1447: Not collective
1449: Input Parameters:
1450: + fvm - The PetscFV object
1451: - sp - The PetscDualSpace object
1453: Level: intermediate
1455: Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1457: .seealso: `PetscFVCreate()`
1458: @*/
1459: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1460: {
1463: PetscDualSpaceDestroy(&fvm->dualSpace);
1464: fvm->dualSpace = sp;
1465: PetscObjectReference((PetscObject)fvm->dualSpace);
1466: return 0;
1467: }
1469: /*@C
1470: PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1472: Not collective
1474: Input Parameter:
1475: . fvm - The PetscFV object
1477: Output Parameter:
1478: . T - The basis function values and derivatives at quadrature points
1480: Note:
1481: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1482: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1483: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1485: Level: intermediate
1487: .seealso: `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1488: @*/
1489: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1490: {
1491: PetscInt npoints;
1492: const PetscReal *points;
1496: PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1497: if (!fvm->T) PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);
1498: *T = fvm->T;
1499: return 0;
1500: }
1502: /*@C
1503: PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1505: Not collective
1507: Input Parameters:
1508: + fvm - The PetscFV object
1509: . nrepl - The number of replicas
1510: . npoints - The number of tabulation points in a replica
1511: . points - The tabulation point coordinates
1512: - K - The order of derivative to tabulate
1514: Output Parameter:
1515: . T - The basis function values and derivative at tabulation points
1517: Note:
1518: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1519: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1520: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e
1522: Level: intermediate
1524: .seealso: `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1525: @*/
1526: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1527: {
1528: PetscInt pdim = 1; /* Dimension of approximation space P */
1529: PetscInt cdim; /* Spatial dimension */
1530: PetscInt Nc; /* Field components */
1531: PetscInt k, p, d, c, e;
1533: if (!npoints || K < 0) {
1534: *T = NULL;
1535: return 0;
1536: }
1540: PetscFVGetSpatialDimension(fvm, &cdim);
1541: PetscFVGetNumComponents(fvm, &Nc);
1542: PetscMalloc1(1, T);
1543: (*T)->K = !cdim ? 0 : K;
1544: (*T)->Nr = nrepl;
1545: (*T)->Np = npoints;
1546: (*T)->Nb = pdim;
1547: (*T)->Nc = Nc;
1548: (*T)->cdim = cdim;
1549: PetscMalloc1((*T)->K + 1, &(*T)->T);
1550: for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
1551: if (K >= 0) {
1552: for (p = 0; p < nrepl * npoints; ++p)
1553: for (d = 0; d < pdim; ++d)
1554: for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1555: }
1556: if (K >= 1) {
1557: for (p = 0; p < nrepl * npoints; ++p)
1558: for (d = 0; d < pdim; ++d)
1559: for (c = 0; c < Nc; ++c)
1560: for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1561: }
1562: if (K >= 2) {
1563: for (p = 0; p < nrepl * npoints; ++p)
1564: for (d = 0; d < pdim; ++d)
1565: for (c = 0; c < Nc; ++c)
1566: for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1567: }
1568: return 0;
1569: }
1571: /*@C
1572: PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1574: Input Parameters:
1575: + fvm - The PetscFV object
1576: . numFaces - The number of cell faces which are not constrained
1577: - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1579: Level: advanced
1581: .seealso: `PetscFVCreate()`
1582: @*/
1583: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1584: {
1586: PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1587: return 0;
1588: }
1590: /*@C
1591: PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1593: Not collective
1595: Input Parameters:
1596: + fvm - The PetscFV object for the field being integrated
1597: . prob - The PetscDS specifing the discretizations and continuum functions
1598: . field - The field being integrated
1599: . Nf - The number of faces in the chunk
1600: . fgeom - The face geometry for each face in the chunk
1601: . neighborVol - The volume for each pair of cells in the chunk
1602: . uL - The state from the cell on the left
1603: - uR - The state from the cell on the right
1605: Output Parameters:
1606: + fluxL - the left fluxes for each face
1607: - fluxR - the right fluxes for each face
1609: Level: developer
1611: .seealso: `PetscFVCreate()`
1612: @*/
1613: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1614: {
1616: PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1617: return 0;
1618: }
1620: /*@
1621: PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1622: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1623: sparsity). It is also used to create an interpolation between regularly refined meshes.
1625: Input Parameter:
1626: . fv - The initial PetscFV
1628: Output Parameter:
1629: . fvRef - The refined PetscFV
1631: Level: advanced
1633: .seealso: `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1634: @*/
1635: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1636: {
1637: PetscDualSpace Q, Qref;
1638: DM K, Kref;
1639: PetscQuadrature q, qref;
1640: DMPolytopeType ct;
1641: DMPlexTransform tr;
1642: PetscReal *v0;
1643: PetscReal *jac, *invjac;
1644: PetscInt numComp, numSubelements, s;
1646: PetscFVGetDualSpace(fv, &Q);
1647: PetscFVGetQuadrature(fv, &q);
1648: PetscDualSpaceGetDM(Q, &K);
1649: /* Create dual space */
1650: PetscDualSpaceDuplicate(Q, &Qref);
1651: DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref);
1652: PetscDualSpaceSetDM(Qref, Kref);
1653: DMDestroy(&Kref);
1654: PetscDualSpaceSetUp(Qref);
1655: /* Create volume */
1656: PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef);
1657: PetscFVSetDualSpace(*fvRef, Qref);
1658: PetscFVGetNumComponents(fv, &numComp);
1659: PetscFVSetNumComponents(*fvRef, numComp);
1660: PetscFVSetUp(*fvRef);
1661: /* Create quadrature */
1662: DMPlexGetCellType(K, 0, &ct);
1663: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1664: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1665: DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1666: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1667: PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1668: for (s = 0; s < numSubelements; ++s) {
1669: PetscQuadrature qs;
1670: const PetscReal *points, *weights;
1671: PetscReal *p, *w;
1672: PetscInt dim, Nc, npoints, np;
1674: PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1675: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1676: np = npoints / numSubelements;
1677: PetscMalloc1(np * dim, &p);
1678: PetscMalloc1(np * Nc, &w);
1679: PetscArraycpy(p, &points[s * np * dim], np * dim);
1680: PetscArraycpy(w, &weights[s * np * Nc], np * Nc);
1681: PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1682: PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1683: PetscQuadratureDestroy(&qs);
1684: }
1685: PetscFVSetQuadrature(*fvRef, qref);
1686: DMPlexTransformDestroy(&tr);
1687: PetscQuadratureDestroy(&qref);
1688: PetscDualSpaceDestroy(&Qref);
1689: return 0;
1690: }
1692: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1693: {
1694: PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;
1696: PetscFree(b);
1697: return 0;
1698: }
1700: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1701: {
1702: PetscInt Nc, c;
1703: PetscViewerFormat format;
1705: PetscFVGetNumComponents(fv, &Nc);
1706: PetscViewerGetFormat(viewer, &format);
1707: PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1708: PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc);
1709: for (c = 0; c < Nc; c++) {
1710: if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1711: }
1712: return 0;
1713: }
1715: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1716: {
1717: PetscBool iascii;
1721: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1722: if (iascii) PetscFVView_Upwind_Ascii(fv, viewer);
1723: return 0;
1724: }
1726: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1727: {
1728: return 0;
1729: }
1731: /*
1732: neighborVol[f*2+0] contains the left geom
1733: neighborVol[f*2+1] contains the right geom
1734: */
1735: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1736: {
1737: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1738: void *rctx;
1739: PetscScalar *flux = fvm->fluxWork;
1740: const PetscScalar *constants;
1741: PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1743: PetscDSGetTotalComponents(prob, &Nc);
1744: PetscDSGetTotalDimension(prob, &totDim);
1745: PetscDSGetFieldOffset(prob, field, &off);
1746: PetscDSGetRiemannSolver(prob, field, &riemann);
1747: PetscDSGetContext(prob, field, &rctx);
1748: PetscDSGetConstants(prob, &numConstants, &constants);
1749: PetscFVGetSpatialDimension(fvm, &dim);
1750: PetscFVGetNumComponents(fvm, &pdim);
1751: for (f = 0; f < Nf; ++f) {
1752: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1753: for (d = 0; d < pdim; ++d) {
1754: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1755: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1756: }
1757: }
1758: return 0;
1759: }
1761: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1762: {
1763: fvm->ops->setfromoptions = NULL;
1764: fvm->ops->setup = PetscFVSetUp_Upwind;
1765: fvm->ops->view = PetscFVView_Upwind;
1766: fvm->ops->destroy = PetscFVDestroy_Upwind;
1767: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1768: return 0;
1769: }
1771: /*MC
1772: PETSCFVUPWIND = "upwind" - A PetscFV object
1774: Level: intermediate
1776: .seealso: `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1777: M*/
1779: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1780: {
1781: PetscFV_Upwind *b;
1784: PetscNew(&b);
1785: fvm->data = b;
1787: PetscFVInitialize_Upwind(fvm);
1788: return 0;
1789: }
1791: #include <petscblaslapack.h>
1793: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1794: {
1795: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
1797: PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1798: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1799: PetscFree(ls);
1800: return 0;
1801: }
1803: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1804: {
1805: PetscInt Nc, c;
1806: PetscViewerFormat format;
1808: PetscFVGetNumComponents(fv, &Nc);
1809: PetscViewerGetFormat(viewer, &format);
1810: PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1811: PetscViewerASCIIPrintf(viewer, " num components: %" PetscInt_FMT "\n", Nc);
1812: for (c = 0; c < Nc; c++) {
1813: if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, " component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1814: }
1815: return 0;
1816: }
1818: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1819: {
1820: PetscBool iascii;
1824: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1825: if (iascii) PetscFVView_LeastSquares_Ascii(fv, viewer);
1826: return 0;
1827: }
1829: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1830: {
1831: return 0;
1832: }
1834: /* Overwrites A. Can only handle full-rank problems with m>=n */
1835: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1836: {
1837: PetscBool debug = PETSC_FALSE;
1838: PetscBLASInt M, N, K, lda, ldb, ldwork, info;
1839: PetscScalar *R, *Q, *Aback, Alpha;
1841: if (debug) {
1842: PetscMalloc1(m * n, &Aback);
1843: PetscArraycpy(Aback, A, m * n);
1844: }
1846: PetscBLASIntCast(m, &M);
1847: PetscBLASIntCast(n, &N);
1848: PetscBLASIntCast(mstride, &lda);
1849: PetscBLASIntCast(worksize, &ldwork);
1850: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1851: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1852: PetscFPTrapPop();
1854: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1856: /* Extract an explicit representation of Q */
1857: Q = Ainv;
1858: PetscArraycpy(Q, A, mstride * n);
1859: K = N; /* full rank */
1860: PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
1863: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1864: Alpha = 1.0;
1865: ldb = lda;
1866: BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1867: /* Ainv is Q, overwritten with inverse */
1869: if (debug) { /* Check that pseudo-inverse worked */
1870: PetscScalar Beta = 0.0;
1871: PetscBLASInt ldc;
1872: K = N;
1873: ldc = N;
1874: BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1875: PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF);
1876: PetscFree(Aback);
1877: }
1878: return 0;
1879: }
1881: /* Overwrites A. Can handle degenerate problems and m<n. */
1882: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1883: {
1884: PetscScalar *Brhs;
1885: PetscScalar *tmpwork;
1886: PetscReal rcond;
1887: #if defined(PETSC_USE_COMPLEX)
1888: PetscInt rworkSize;
1889: PetscReal *rwork;
1890: #endif
1891: PetscInt i, j, maxmn;
1892: PetscBLASInt M, N, lda, ldb, ldwork;
1893: PetscBLASInt nrhs, irank, info;
1895: /* initialize to identity */
1896: tmpwork = work;
1897: Brhs = Ainv;
1898: maxmn = PetscMax(m, n);
1899: for (j = 0; j < maxmn; j++) {
1900: for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
1901: }
1903: PetscBLASIntCast(m, &M);
1904: PetscBLASIntCast(n, &N);
1905: PetscBLASIntCast(mstride, &lda);
1906: PetscBLASIntCast(maxmn, &ldb);
1907: PetscBLASIntCast(worksize, &ldwork);
1908: rcond = -1;
1909: nrhs = M;
1910: #if defined(PETSC_USE_COMPLEX)
1911: rworkSize = 5 * PetscMin(M, N);
1912: PetscMalloc1(rworkSize, &rwork);
1913: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1914: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
1915: PetscFPTrapPop();
1916: PetscFree(rwork);
1917: #else
1918: nrhs = M;
1919: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1920: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
1921: PetscFPTrapPop();
1922: #endif
1924: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1926: return 0;
1927: }
1929: #if 0
1930: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1931: {
1932: PetscReal grad[2] = {0, 0};
1933: const PetscInt *faces;
1934: PetscInt numFaces, f;
1936: DMPlexGetConeSize(dm, cell, &numFaces);
1937: DMPlexGetCone(dm, cell, &faces);
1938: for (f = 0; f < numFaces; ++f) {
1939: const PetscInt *fcells;
1940: const CellGeom *cg1;
1941: const FaceGeom *fg;
1943: DMPlexGetSupport(dm, faces[f], &fcells);
1944: DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1945: for (i = 0; i < 2; ++i) {
1946: PetscScalar du;
1948: if (fcells[i] == c) continue;
1949: DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1950: du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1951: grad[0] += fg->grad[!i][0] * du;
1952: grad[1] += fg->grad[!i][1] * du;
1953: }
1954: }
1955: PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1956: return 0;
1957: }
1958: #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: static 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;
1979: if (numFaces > maxFaces) {
1981: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
1982: }
1983: PetscFVGetSpatialDimension(fvm, &dim);
1984: for (f = 0; f < numFaces; ++f) {
1985: for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
1986: }
1987: /* Overwrites B with garbage, returns Binv in row-major format */
1988: if (useSVD) {
1989: PetscInt maxmn = PetscMax(numFaces, dim);
1990: PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1991: /* Binv shaped in column-major, coldim=maxmn.*/
1992: for (f = 0; f < numFaces; ++f) {
1993: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
1994: }
1995: } else {
1996: PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1997: /* Binv shaped in row-major, rowdim=maxFaces.*/
1998: for (f = 0; f < numFaces; ++f) {
1999: for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2000: }
2001: }
2002: return 0;
2003: }
2005: /*
2006: neighborVol[f*2+0] contains the left geom
2007: neighborVol[f*2+1] contains the right geom
2008: */
2009: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2010: {
2011: void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2012: void *rctx;
2013: PetscScalar *flux = fvm->fluxWork;
2014: const PetscScalar *constants;
2015: PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2017: PetscDSGetTotalComponents(prob, &Nc);
2018: PetscDSGetTotalDimension(prob, &totDim);
2019: PetscDSGetFieldOffset(prob, field, &off);
2020: PetscDSGetRiemannSolver(prob, field, &riemann);
2021: PetscDSGetContext(prob, field, &rctx);
2022: PetscDSGetConstants(prob, &numConstants, &constants);
2023: PetscFVGetSpatialDimension(fvm, &dim);
2024: PetscFVGetNumComponents(fvm, &pdim);
2025: for (f = 0; f < Nf; ++f) {
2026: (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2027: for (d = 0; d < pdim; ++d) {
2028: fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2029: fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2030: }
2031: }
2032: return 0;
2033: }
2035: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2036: {
2037: PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2038: PetscInt dim, m, n, nrhs, minmn, maxmn;
2041: PetscFVGetSpatialDimension(fvm, &dim);
2042: PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2043: ls->maxFaces = maxFaces;
2044: m = ls->maxFaces;
2045: n = dim;
2046: nrhs = ls->maxFaces;
2047: minmn = PetscMin(m, n);
2048: maxmn = PetscMax(m, n);
2049: ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2050: PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work);
2051: return 0;
2052: }
2054: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2055: {
2056: fvm->ops->setfromoptions = NULL;
2057: fvm->ops->setup = PetscFVSetUp_LeastSquares;
2058: fvm->ops->view = PetscFVView_LeastSquares;
2059: fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2060: fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2061: fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2062: return 0;
2063: }
2065: /*MC
2066: PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2068: Level: intermediate
2070: .seealso: `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2071: M*/
2073: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2074: {
2075: PetscFV_LeastSquares *ls;
2078: PetscNew(&ls);
2079: fvm->data = ls;
2081: ls->maxFaces = -1;
2082: ls->workSize = -1;
2083: ls->B = NULL;
2084: ls->Binv = NULL;
2085: ls->tau = NULL;
2086: ls->work = NULL;
2088: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2089: PetscFVInitialize_LeastSquares(fvm);
2090: PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2091: return 0;
2092: }
2094: /*@
2095: PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2097: Not collective
2099: Input parameters:
2100: + fvm - The PetscFV object
2101: - maxFaces - The maximum number of cell faces
2103: Level: intermediate
2105: .seealso: `PetscFVCreate()`, `PETSCFVLEASTSQUARES`
2106: @*/
2107: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2108: {
2110: PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2111: return 0;
2112: }