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: }