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: + sname    - The name of a new user-defined creation routine
 26: - function - The creation routine

 28:   Example Usage:
 29: .vb
 30:     PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
 31: .ve

 33:   Then, your `PetscLimiter` type can be chosen with the procedural interface via
 34: .vb
 35:     PetscLimiterCreate(MPI_Comm, PetscLimiter *);
 36:     PetscLimiterSetType(PetscLimiter, "my_lim");
 37: .ve
 38:   or at runtime via the option
 39: .vb
 40:     -petsclimiter_type my_lim
 41: .ve

 43:   Level: advanced

 45:   Note:
 46:   `PetscLimiterRegister()` may be called multiple times to add several user-defined PetscLimiters

 48: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterRegisterAll()`, `PetscLimiterRegisterDestroy()`
 49: @*/
 50: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
 51: {
 52:   PetscFunctionBegin;
 53:   PetscCall(PetscFunctionListAdd(&PetscLimiterList, sname, function));
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 55: }

 57: /*@C
 58:   PetscLimiterSetType - Builds a `PetscLimiter` for a given `PetscLimiterType`

 60:   Collective

 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: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
 72: @*/
 73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
 74: {
 75:   PetscErrorCode (*r)(PetscLimiter);
 76:   PetscBool match;

 78:   PetscFunctionBegin;
 80:   PetscCall(PetscObjectTypeCompare((PetscObject)lim, name, &match));
 81:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

 83:   PetscCall(PetscLimiterRegisterAll());
 84:   PetscCall(PetscFunctionListFind(PetscLimiterList, name, &r));
 85:   PetscCheck(r, PetscObjectComm((PetscObject)lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);

 87:   PetscTryTypeMethod(lim, destroy);
 88:   lim->ops->destroy = NULL;

 90:   PetscCall((*r)(lim));
 91:   PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: /*@C
 96:   PetscLimiterGetType - Gets the `PetscLimiterType` name (as a string) from the `PetscLimiter`.

 98:   Not Collective

100:   Input Parameter:
101: . lim - The `PetscLimiter`

103:   Output Parameter:
104: . name - The `PetscLimiterType`

106:   Level: intermediate

108: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
112:   PetscFunctionBegin;
114:   PetscAssertPointer(name, 2);
115:   PetscCall(PetscLimiterRegisterAll());
116:   *name = ((PetscObject)lim)->type_name;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /*@C
121:   PetscLimiterViewFromOptions - View a `PetscLimiter` based on values in the options database

123:   Collective

125:   Input Parameters:
126: + A    - the `PetscLimiter` object to view
127: . obj  - Optional object that provides the options prefix to use
128: - name - command line option name

130:   Level: intermediate

132: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
133: @*/
134: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
135: {
136:   PetscFunctionBegin;
138:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: /*@C
143:   PetscLimiterView - Views a `PetscLimiter`

145:   Collective

147:   Input Parameters:
148: + lim - the `PetscLimiter` object to view
149: - v   - the viewer

151:   Level: beginner

153: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
154: @*/
155: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
156: {
157:   PetscFunctionBegin;
159:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v));
160:   PetscTryTypeMethod(lim, view, v);
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: /*@
165:   PetscLimiterSetFromOptions - sets parameters in a `PetscLimiter` from the options database

167:   Collective

169:   Input Parameter:
170: . lim - the `PetscLimiter` object to set options for

172:   Level: intermediate

174: .seealso: `PetscLimiter`, `PetscLimiterView()`
175: @*/
176: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
177: {
178:   const char *defaultType;
179:   char        name[256];
180:   PetscBool   flg;

182:   PetscFunctionBegin;
184:   if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
185:   else defaultType = ((PetscObject)lim)->type_name;
186:   PetscCall(PetscLimiterRegisterAll());

188:   PetscObjectOptionsBegin((PetscObject)lim);
189:   PetscCall(PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg));
190:   if (flg) {
191:     PetscCall(PetscLimiterSetType(lim, name));
192:   } else if (!((PetscObject)lim)->type_name) {
193:     PetscCall(PetscLimiterSetType(lim, defaultType));
194:   }
195:   PetscTryTypeMethod(lim, setfromoptions);
196:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
197:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject));
198:   PetscOptionsEnd();
199:   PetscCall(PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view"));
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }

203: /*@C
204:   PetscLimiterSetUp - Construct data structures for the `PetscLimiter`

206:   Collective

208:   Input Parameter:
209: . lim - the `PetscLimiter` object to setup

211:   Level: intermediate

213: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscLimiterDestroy()`
214: @*/
215: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
216: {
217:   PetscFunctionBegin;
219:   PetscTryTypeMethod(lim, setup);
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: /*@
224:   PetscLimiterDestroy - Destroys a `PetscLimiter` object

226:   Collective

228:   Input Parameter:
229: . lim - the `PetscLimiter` object to destroy

231:   Level: beginner

233: .seealso: `PetscLimiter`, `PetscLimiterView()`
234: @*/
235: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
236: {
237:   PetscFunctionBegin;
238:   if (!*lim) PetscFunctionReturn(PETSC_SUCCESS);

241:   if (--((PetscObject)*lim)->refct > 0) {
242:     *lim = NULL;
243:     PetscFunctionReturn(PETSC_SUCCESS);
244:   }
245:   ((PetscObject)*lim)->refct = 0;

247:   PetscTryTypeMethod(*lim, destroy);
248:   PetscCall(PetscHeaderDestroy(lim));
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: /*@
253:   PetscLimiterCreate - Creates an empty `PetscLimiter` object. The type can then be set with `PetscLimiterSetType()`.

255:   Collective

257:   Input Parameter:
258: . comm - The communicator for the `PetscLimiter` object

260:   Output Parameter:
261: . lim - The `PetscLimiter` object

263:   Level: beginner

265: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
266: @*/
267: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
268: {
269:   PetscLimiter l;

271:   PetscFunctionBegin;
272:   PetscAssertPointer(lim, 2);
273:   PetscCall(PetscCitationsRegister(LimiterCitation, &Limitercite));
274:   *lim = NULL;
275:   PetscCall(PetscFVInitializePackage());

277:   PetscCall(PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView));

279:   *lim = l;
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: /*@
284:   PetscLimiterLimit - Limit the flux

286:   Input Parameters:
287: + lim  - The `PetscLimiter`
288: - flim - The input field

290:   Output Parameter:
291: . phi - The limited field

293:   Level: beginner

295:   Note:
296:   Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
297: .vb
298:  The classical flux-limited formulation is psi(r) where

300:  r = (u[0] - u[-1]) / (u[1] - u[0])

302:  The second order TVD region is bounded by

304:  psi_minmod(r) = min(r,1)      and        psi_superbee(r) = min(2, 2r, max(1,r))

306:  where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
307:  phi(r)(r+1)/2 in which the reconstructed interface values are

309:  u(v) = u[0] + phi(r) (grad u)[0] v

311:  where v is the vector from centroid to quadrature point. In these variables, the usual limiters become

313:  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))

315:  For a nicer symmetric formulation, rewrite in terms of

317:  f = (u[0] - u[-1]) / (u[1] - u[-1])

319:  where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition

321:  phi(r) = phi(1/r)

323:  becomes

325:  w(f) = w(1-f).

327:  The limiters below implement this final form w(f). The reference methods are

329:  w_minmod(f) = 2 min(f,(1-f))             w_superbee(r) = 4 min((1-f), f)
330: .ve

332: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
333: @*/
334: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
335: {
336:   PetscFunctionBegin;
338:   PetscAssertPointer(phi, 3);
339:   PetscUseTypeMethod(lim, limit, flim, phi);
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
344: {
345:   PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;

347:   PetscFunctionBegin;
348:   PetscCall(PetscFree(l));
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

352: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
353: {
354:   PetscViewerFormat format;

356:   PetscFunctionBegin;
357:   PetscCall(PetscViewerGetFormat(viewer, &format));
358:   PetscCall(PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n"));
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
363: {
364:   PetscBool iascii;

366:   PetscFunctionBegin;
369:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
370:   if (iascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
375: {
376:   PetscFunctionBegin;
377:   *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
382: {
383:   PetscFunctionBegin;
384:   lim->ops->view    = PetscLimiterView_Sin;
385:   lim->ops->destroy = PetscLimiterDestroy_Sin;
386:   lim->ops->limit   = PetscLimiterLimit_Sin;
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: /*MC
391:   PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation

393:   Level: intermediate

395: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
396: M*/

398: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
399: {
400:   PetscLimiter_Sin *l;

402:   PetscFunctionBegin;
404:   PetscCall(PetscNew(&l));
405:   lim->data = l;

407:   PetscCall(PetscLimiterInitialize_Sin(lim));
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }

411: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
412: {
413:   PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;

415:   PetscFunctionBegin;
416:   PetscCall(PetscFree(l));
417:   PetscFunctionReturn(PETSC_SUCCESS);
418: }

420: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
421: {
422:   PetscViewerFormat format;

424:   PetscFunctionBegin;
425:   PetscCall(PetscViewerGetFormat(viewer, &format));
426:   PetscCall(PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n"));
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
431: {
432:   PetscBool iascii;

434:   PetscFunctionBegin;
437:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
438:   if (iascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
443: {
444:   PetscFunctionBegin;
445:   *phi = 0.0;
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
450: {
451:   PetscFunctionBegin;
452:   lim->ops->view    = PetscLimiterView_Zero;
453:   lim->ops->destroy = PetscLimiterDestroy_Zero;
454:   lim->ops->limit   = PetscLimiterLimit_Zero;
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: /*MC
459:   PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation

461:   Level: intermediate

463: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
464: M*/

466: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
467: {
468:   PetscLimiter_Zero *l;

470:   PetscFunctionBegin;
472:   PetscCall(PetscNew(&l));
473:   lim->data = l;

475:   PetscCall(PetscLimiterInitialize_Zero(lim));
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
480: {
481:   PetscLimiter_None *l = (PetscLimiter_None *)lim->data;

483:   PetscFunctionBegin;
484:   PetscCall(PetscFree(l));
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
489: {
490:   PetscViewerFormat format;

492:   PetscFunctionBegin;
493:   PetscCall(PetscViewerGetFormat(viewer, &format));
494:   PetscCall(PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n"));
495:   PetscFunctionReturn(PETSC_SUCCESS);
496: }

498: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
499: {
500:   PetscBool iascii;

502:   PetscFunctionBegin;
505:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
506:   if (iascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

510: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
511: {
512:   PetscFunctionBegin;
513:   *phi = 1.0;
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
518: {
519:   PetscFunctionBegin;
520:   lim->ops->view    = PetscLimiterView_None;
521:   lim->ops->destroy = PetscLimiterDestroy_None;
522:   lim->ops->limit   = PetscLimiterLimit_None;
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: /*MC
527:   PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation

529:   Level: intermediate

531: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
532: M*/

534: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
535: {
536:   PetscLimiter_None *l;

538:   PetscFunctionBegin;
540:   PetscCall(PetscNew(&l));
541:   lim->data = l;

543:   PetscCall(PetscLimiterInitialize_None(lim));
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
548: {
549:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;

551:   PetscFunctionBegin;
552:   PetscCall(PetscFree(l));
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
557: {
558:   PetscViewerFormat format;

560:   PetscFunctionBegin;
561:   PetscCall(PetscViewerGetFormat(viewer, &format));
562:   PetscCall(PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n"));
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
567: {
568:   PetscBool iascii;

570:   PetscFunctionBegin;
573:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
574:   if (iascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
575:   PetscFunctionReturn(PETSC_SUCCESS);
576: }

578: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
579: {
580:   PetscFunctionBegin;
581:   *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
582:   PetscFunctionReturn(PETSC_SUCCESS);
583: }

585: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
586: {
587:   PetscFunctionBegin;
588:   lim->ops->view    = PetscLimiterView_Minmod;
589:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
590:   lim->ops->limit   = PetscLimiterLimit_Minmod;
591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }

594: /*MC
595:   PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation

597:   Level: intermediate

599: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
600: M*/

602: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
603: {
604:   PetscLimiter_Minmod *l;

606:   PetscFunctionBegin;
608:   PetscCall(PetscNew(&l));
609:   lim->data = l;

611:   PetscCall(PetscLimiterInitialize_Minmod(lim));
612:   PetscFunctionReturn(PETSC_SUCCESS);
613: }

615: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
616: {
617:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;

619:   PetscFunctionBegin;
620:   PetscCall(PetscFree(l));
621:   PetscFunctionReturn(PETSC_SUCCESS);
622: }

624: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
625: {
626:   PetscViewerFormat format;

628:   PetscFunctionBegin;
629:   PetscCall(PetscViewerGetFormat(viewer, &format));
630:   PetscCall(PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n"));
631:   PetscFunctionReturn(PETSC_SUCCESS);
632: }

634: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
635: {
636:   PetscBool iascii;

638:   PetscFunctionBegin;
641:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
642:   if (iascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
643:   PetscFunctionReturn(PETSC_SUCCESS);
644: }

646: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
647: {
648:   PetscFunctionBegin;
649:   *phi = PetscMax(0, 4 * f * (1 - f));
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
654: {
655:   PetscFunctionBegin;
656:   lim->ops->view    = PetscLimiterView_VanLeer;
657:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
658:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
659:   PetscFunctionReturn(PETSC_SUCCESS);
660: }

662: /*MC
663:   PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation

665:   Level: intermediate

667: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
668: M*/

670: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
671: {
672:   PetscLimiter_VanLeer *l;

674:   PetscFunctionBegin;
676:   PetscCall(PetscNew(&l));
677:   lim->data = l;

679:   PetscCall(PetscLimiterInitialize_VanLeer(lim));
680:   PetscFunctionReturn(PETSC_SUCCESS);
681: }

683: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
684: {
685:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;

687:   PetscFunctionBegin;
688:   PetscCall(PetscFree(l));
689:   PetscFunctionReturn(PETSC_SUCCESS);
690: }

692: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
693: {
694:   PetscViewerFormat format;

696:   PetscFunctionBegin;
697:   PetscCall(PetscViewerGetFormat(viewer, &format));
698:   PetscCall(PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n"));
699:   PetscFunctionReturn(PETSC_SUCCESS);
700: }

702: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
703: {
704:   PetscBool iascii;

706:   PetscFunctionBegin;
709:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
710:   if (iascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
711:   PetscFunctionReturn(PETSC_SUCCESS);
712: }

714: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
715: {
716:   PetscFunctionBegin;
717:   *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
718:   PetscFunctionReturn(PETSC_SUCCESS);
719: }

721: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
722: {
723:   PetscFunctionBegin;
724:   lim->ops->view    = PetscLimiterView_VanAlbada;
725:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
726:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
727:   PetscFunctionReturn(PETSC_SUCCESS);
728: }

730: /*MC
731:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation

733:   Level: intermediate

735: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
736: M*/

738: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
739: {
740:   PetscLimiter_VanAlbada *l;

742:   PetscFunctionBegin;
744:   PetscCall(PetscNew(&l));
745:   lim->data = l;

747:   PetscCall(PetscLimiterInitialize_VanAlbada(lim));
748:   PetscFunctionReturn(PETSC_SUCCESS);
749: }

751: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
752: {
753:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;

755:   PetscFunctionBegin;
756:   PetscCall(PetscFree(l));
757:   PetscFunctionReturn(PETSC_SUCCESS);
758: }

760: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
761: {
762:   PetscViewerFormat format;

764:   PetscFunctionBegin;
765:   PetscCall(PetscViewerGetFormat(viewer, &format));
766:   PetscCall(PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n"));
767:   PetscFunctionReturn(PETSC_SUCCESS);
768: }

770: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
771: {
772:   PetscBool iascii;

774:   PetscFunctionBegin;
777:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
778:   if (iascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
779:   PetscFunctionReturn(PETSC_SUCCESS);
780: }

782: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
783: {
784:   PetscFunctionBegin;
785:   *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
786:   PetscFunctionReturn(PETSC_SUCCESS);
787: }

789: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
790: {
791:   PetscFunctionBegin;
792:   lim->ops->view    = PetscLimiterView_Superbee;
793:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
794:   lim->ops->limit   = PetscLimiterLimit_Superbee;
795:   PetscFunctionReturn(PETSC_SUCCESS);
796: }

798: /*MC
799:   PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation

801:   Level: intermediate

803: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
804: M*/

806: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
807: {
808:   PetscLimiter_Superbee *l;

810:   PetscFunctionBegin;
812:   PetscCall(PetscNew(&l));
813:   lim->data = l;

815:   PetscCall(PetscLimiterInitialize_Superbee(lim));
816:   PetscFunctionReturn(PETSC_SUCCESS);
817: }

819: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
820: {
821:   PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;

823:   PetscFunctionBegin;
824:   PetscCall(PetscFree(l));
825:   PetscFunctionReturn(PETSC_SUCCESS);
826: }

828: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
829: {
830:   PetscViewerFormat format;

832:   PetscFunctionBegin;
833:   PetscCall(PetscViewerGetFormat(viewer, &format));
834:   PetscCall(PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n"));
835:   PetscFunctionReturn(PETSC_SUCCESS);
836: }

838: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
839: {
840:   PetscBool iascii;

842:   PetscFunctionBegin;
845:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
846:   if (iascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
847:   PetscFunctionReturn(PETSC_SUCCESS);
848: }

850: /* aka Barth-Jespersen */
851: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
852: {
853:   PetscFunctionBegin;
854:   *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
855:   PetscFunctionReturn(PETSC_SUCCESS);
856: }

858: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
859: {
860:   PetscFunctionBegin;
861:   lim->ops->view    = PetscLimiterView_MC;
862:   lim->ops->destroy = PetscLimiterDestroy_MC;
863:   lim->ops->limit   = PetscLimiterLimit_MC;
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }

867: /*MC
868:   PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation

870:   Level: intermediate

872: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
873: M*/

875: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
876: {
877:   PetscLimiter_MC *l;

879:   PetscFunctionBegin;
881:   PetscCall(PetscNew(&l));
882:   lim->data = l;

884:   PetscCall(PetscLimiterInitialize_MC(lim));
885:   PetscFunctionReturn(PETSC_SUCCESS);
886: }

888: PetscClassId PETSCFV_CLASSID = 0;

890: PetscFunctionList PetscFVList              = NULL;
891: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

893: /*@C
894:   PetscFVRegister - Adds a new `PetscFV` implementation

896:   Not Collective

898:   Input Parameters:
899: + sname    - The name of a new user-defined creation routine
900: - function - The creation routine itself

902:   Example Usage:
903: .vb
904:     PetscFVRegister("my_fv", MyPetscFVCreate);
905: .ve

907:   Then, your PetscFV type can be chosen with the procedural interface via
908: .vb
909:     PetscFVCreate(MPI_Comm, PetscFV *);
910:     PetscFVSetType(PetscFV, "my_fv");
911: .ve
912:   or at runtime via the option
913: .vb
914:     -petscfv_type my_fv
915: .ve

917:   Level: advanced

919:   Note:
920:   `PetscFVRegister()` may be called multiple times to add several user-defined PetscFVs

922: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
923: @*/
924: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
925: {
926:   PetscFunctionBegin;
927:   PetscCall(PetscFunctionListAdd(&PetscFVList, sname, function));
928:   PetscFunctionReturn(PETSC_SUCCESS);
929: }

931: /*@C
932:   PetscFVSetType - Builds a particular `PetscFV`

934:   Collective

936:   Input Parameters:
937: + fvm  - The `PetscFV` object
938: - name - The type of FVM space

940:   Options Database Key:
941: . -petscfv_type <type> - Sets the `PetscFVType`; use -help for a list of available types

943:   Level: intermediate

945: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
946: @*/
947: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
948: {
949:   PetscErrorCode (*r)(PetscFV);
950:   PetscBool match;

952:   PetscFunctionBegin;
954:   PetscCall(PetscObjectTypeCompare((PetscObject)fvm, name, &match));
955:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

957:   PetscCall(PetscFVRegisterAll());
958:   PetscCall(PetscFunctionListFind(PetscFVList, name, &r));
959:   PetscCheck(r, PetscObjectComm((PetscObject)fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);

961:   PetscTryTypeMethod(fvm, destroy);
962:   fvm->ops->destroy = NULL;

964:   PetscCall((*r)(fvm));
965:   PetscCall(PetscObjectChangeTypeName((PetscObject)fvm, name));
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: /*@C
970:   PetscFVGetType - Gets the `PetscFVType` (as a string) from a `PetscFV`.

972:   Not Collective

974:   Input Parameter:
975: . fvm - The `PetscFV`

977:   Output Parameter:
978: . name - The `PetscFVType` name

980:   Level: intermediate

982: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
983: @*/
984: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
985: {
986:   PetscFunctionBegin;
988:   PetscAssertPointer(name, 2);
989:   PetscCall(PetscFVRegisterAll());
990:   *name = ((PetscObject)fvm)->type_name;
991:   PetscFunctionReturn(PETSC_SUCCESS);
992: }

994: /*@C
995:   PetscFVViewFromOptions - View a `PetscFV` based on values in the options database

997:   Collective

999:   Input Parameters:
1000: + A    - the `PetscFV` object
1001: . obj  - Optional object that provides the options prefix
1002: - name - command line option name

1004:   Level: intermediate

1006: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
1007: @*/
1008: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
1009: {
1010:   PetscFunctionBegin;
1012:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
1013:   PetscFunctionReturn(PETSC_SUCCESS);
1014: }

1016: /*@C
1017:   PetscFVView - Views a `PetscFV`

1019:   Collective

1021:   Input Parameters:
1022: + fvm - the `PetscFV` object to view
1023: - v   - the viewer

1025:   Level: beginner

1027: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
1028: @*/
1029: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1030: {
1031:   PetscFunctionBegin;
1033:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v));
1034:   PetscTryTypeMethod(fvm, view, v);
1035:   PetscFunctionReturn(PETSC_SUCCESS);
1036: }

1038: /*@
1039:   PetscFVSetFromOptions - sets parameters in a `PetscFV` from the options database

1041:   Collective

1043:   Input Parameter:
1044: . fvm - the `PetscFV` object to set options for

1046:   Options Database Key:
1047: . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated

1049:   Level: intermediate

1051: .seealso: `PetscFV`, `PetscFVView()`
1052: @*/
1053: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1054: {
1055:   const char *defaultType;
1056:   char        name[256];
1057:   PetscBool   flg;

1059:   PetscFunctionBegin;
1061:   if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
1062:   else defaultType = ((PetscObject)fvm)->type_name;
1063:   PetscCall(PetscFVRegisterAll());

1065:   PetscObjectOptionsBegin((PetscObject)fvm);
1066:   PetscCall(PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg));
1067:   if (flg) {
1068:     PetscCall(PetscFVSetType(fvm, name));
1069:   } else if (!((PetscObject)fvm)->type_name) {
1070:     PetscCall(PetscFVSetType(fvm, defaultType));
1071:   }
1072:   PetscCall(PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL));
1073:   PetscTryTypeMethod(fvm, setfromoptions);
1074:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1075:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject));
1076:   PetscCall(PetscLimiterSetFromOptions(fvm->limiter));
1077:   PetscOptionsEnd();
1078:   PetscCall(PetscFVViewFromOptions(fvm, NULL, "-petscfv_view"));
1079:   PetscFunctionReturn(PETSC_SUCCESS);
1080: }

1082: /*@
1083:   PetscFVSetUp - Setup the data structures for the `PetscFV` based on the `PetscFVType` provided by `PetscFVSetType()`

1085:   Collective

1087:   Input Parameter:
1088: . fvm - the `PetscFV` object to setup

1090:   Level: intermediate

1092: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1093: @*/
1094: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1095: {
1096:   PetscFunctionBegin;
1098:   PetscCall(PetscLimiterSetUp(fvm->limiter));
1099:   PetscTryTypeMethod(fvm, setup);
1100:   PetscFunctionReturn(PETSC_SUCCESS);
1101: }

1103: /*@
1104:   PetscFVDestroy - Destroys a `PetscFV` object

1106:   Collective

1108:   Input Parameter:
1109: . fvm - the `PetscFV` object to destroy

1111:   Level: beginner

1113: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1114: @*/
1115: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1116: {
1117:   PetscInt i;

1119:   PetscFunctionBegin;
1120:   if (!*fvm) PetscFunctionReturn(PETSC_SUCCESS);

1123:   if (--((PetscObject)*fvm)->refct > 0) {
1124:     *fvm = NULL;
1125:     PetscFunctionReturn(PETSC_SUCCESS);
1126:   }
1127:   ((PetscObject)*fvm)->refct = 0;

1129:   for (i = 0; i < (*fvm)->numComponents; i++) PetscCall(PetscFree((*fvm)->componentNames[i]));
1130:   PetscCall(PetscFree((*fvm)->componentNames));
1131:   PetscCall(PetscLimiterDestroy(&(*fvm)->limiter));
1132:   PetscCall(PetscDualSpaceDestroy(&(*fvm)->dualSpace));
1133:   PetscCall(PetscFree((*fvm)->fluxWork));
1134:   PetscCall(PetscQuadratureDestroy(&(*fvm)->quadrature));
1135:   PetscCall(PetscTabulationDestroy(&(*fvm)->T));

1137:   PetscTryTypeMethod(*fvm, destroy);
1138:   PetscCall(PetscHeaderDestroy(fvm));
1139:   PetscFunctionReturn(PETSC_SUCCESS);
1140: }

1142: /*@
1143:   PetscFVCreate - Creates an empty `PetscFV` object. The type can then be set with `PetscFVSetType()`.

1145:   Collective

1147:   Input Parameter:
1148: . comm - The communicator for the `PetscFV` object

1150:   Output Parameter:
1151: . fvm - The `PetscFV` object

1153:   Level: beginner

1155: .seealso: `PetscFVSetUp()`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1156: @*/
1157: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1158: {
1159:   PetscFV f;

1161:   PetscFunctionBegin;
1162:   PetscAssertPointer(fvm, 2);
1163:   *fvm = NULL;
1164:   PetscCall(PetscFVInitializePackage());

1166:   PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
1167:   PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));

1169:   PetscCall(PetscLimiterCreate(comm, &f->limiter));
1170:   f->numComponents    = 1;
1171:   f->dim              = 0;
1172:   f->computeGradients = PETSC_FALSE;
1173:   f->fluxWork         = NULL;
1174:   PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));

1176:   *fvm = f;
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: /*@
1181:   PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`

1183:   Logically Collective

1185:   Input Parameters:
1186: + fvm - the `PetscFV` object
1187: - lim - The `PetscLimiter`

1189:   Level: intermediate

1191: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1192: @*/
1193: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1194: {
1195:   PetscFunctionBegin;
1198:   PetscCall(PetscLimiterDestroy(&fvm->limiter));
1199:   PetscCall(PetscObjectReference((PetscObject)lim));
1200:   fvm->limiter = lim;
1201:   PetscFunctionReturn(PETSC_SUCCESS);
1202: }

1204: /*@
1205:   PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`

1207:   Not Collective

1209:   Input Parameter:
1210: . fvm - the `PetscFV` object

1212:   Output Parameter:
1213: . lim - The `PetscLimiter`

1215:   Level: intermediate

1217: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1218: @*/
1219: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1220: {
1221:   PetscFunctionBegin;
1223:   PetscAssertPointer(lim, 2);
1224:   *lim = fvm->limiter;
1225:   PetscFunctionReturn(PETSC_SUCCESS);
1226: }

1228: /*@
1229:   PetscFVSetNumComponents - Set the number of field components in a `PetscFV`

1231:   Logically Collective

1233:   Input Parameters:
1234: + fvm  - the `PetscFV` object
1235: - comp - The number of components

1237:   Level: intermediate

1239: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1240: @*/
1241: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1242: {
1243:   PetscFunctionBegin;
1245:   if (fvm->numComponents != comp) {
1246:     PetscInt i;

1248:     for (i = 0; i < fvm->numComponents; i++) PetscCall(PetscFree(fvm->componentNames[i]));
1249:     PetscCall(PetscFree(fvm->componentNames));
1250:     PetscCall(PetscCalloc1(comp, &fvm->componentNames));
1251:   }
1252:   fvm->numComponents = comp;
1253:   PetscCall(PetscFree(fvm->fluxWork));
1254:   PetscCall(PetscMalloc1(comp, &fvm->fluxWork));
1255:   PetscFunctionReturn(PETSC_SUCCESS);
1256: }

1258: /*@
1259:   PetscFVGetNumComponents - Get the number of field components in a `PetscFV`

1261:   Not Collective

1263:   Input Parameter:
1264: . fvm - the `PetscFV` object

1266:   Output Parameter:
1267: . comp - The number of components

1269:   Level: intermediate

1271: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1272: @*/
1273: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1274: {
1275:   PetscFunctionBegin;
1277:   PetscAssertPointer(comp, 2);
1278:   *comp = fvm->numComponents;
1279:   PetscFunctionReturn(PETSC_SUCCESS);
1280: }

1282: /*@C
1283:   PetscFVSetComponentName - Set the name of a component (used in output and viewing) in a `PetscFV`

1285:   Logically Collective

1287:   Input Parameters:
1288: + fvm  - the `PetscFV` object
1289: . comp - the component number
1290: - name - the component name

1292:   Level: intermediate

1294: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1295: @*/
1296: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1297: {
1298:   PetscFunctionBegin;
1299:   PetscCall(PetscFree(fvm->componentNames[comp]));
1300:   PetscCall(PetscStrallocpy(name, &fvm->componentNames[comp]));
1301:   PetscFunctionReturn(PETSC_SUCCESS);
1302: }

1304: /*@C
1305:   PetscFVGetComponentName - Get the name of a component (used in output and viewing) in a `PetscFV`

1307:   Logically Collective

1309:   Input Parameters:
1310: + fvm  - the `PetscFV` object
1311: - comp - the component number

1313:   Output Parameter:
1314: . name - the component name

1316:   Level: intermediate

1318: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1319: @*/
1320: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1321: {
1322:   PetscFunctionBegin;
1323:   *name = fvm->componentNames[comp];
1324:   PetscFunctionReturn(PETSC_SUCCESS);
1325: }

1327: /*@
1328:   PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`

1330:   Logically Collective

1332:   Input Parameters:
1333: + fvm - the `PetscFV` object
1334: - dim - The spatial dimension

1336:   Level: intermediate

1338: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
1339: @*/
1340: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1341: {
1342:   PetscFunctionBegin;
1344:   fvm->dim = dim;
1345:   PetscFunctionReturn(PETSC_SUCCESS);
1346: }

1348: /*@
1349:   PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`

1351:   Not Collective

1353:   Input Parameter:
1354: . fvm - the `PetscFV` object

1356:   Output Parameter:
1357: . dim - The spatial dimension

1359:   Level: intermediate

1361: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1362: @*/
1363: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1364: {
1365:   PetscFunctionBegin;
1367:   PetscAssertPointer(dim, 2);
1368:   *dim = fvm->dim;
1369:   PetscFunctionReturn(PETSC_SUCCESS);
1370: }

1372: /*@
1373:   PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`

1375:   Logically Collective

1377:   Input Parameters:
1378: + fvm              - the `PetscFV` object
1379: - computeGradients - Flag to compute cell gradients

1381:   Level: intermediate

1383: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1384: @*/
1385: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1386: {
1387:   PetscFunctionBegin;
1389:   fvm->computeGradients = computeGradients;
1390:   PetscFunctionReturn(PETSC_SUCCESS);
1391: }

1393: /*@
1394:   PetscFVGetComputeGradients - Return flag for computation of cell gradients on a `PetscFV`

1396:   Not Collective

1398:   Input Parameter:
1399: . fvm - the `PetscFV` object

1401:   Output Parameter:
1402: . computeGradients - Flag to compute cell gradients

1404:   Level: intermediate

1406: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1407: @*/
1408: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1409: {
1410:   PetscFunctionBegin;
1412:   PetscAssertPointer(computeGradients, 2);
1413:   *computeGradients = fvm->computeGradients;
1414:   PetscFunctionReturn(PETSC_SUCCESS);
1415: }

1417: /*@
1418:   PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`

1420:   Logically Collective

1422:   Input Parameters:
1423: + fvm - the `PetscFV` object
1424: - q   - The `PetscQuadrature`

1426:   Level: intermediate

1428: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1429: @*/
1430: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1431: {
1432:   PetscFunctionBegin;
1434:   PetscCall(PetscObjectReference((PetscObject)q));
1435:   PetscCall(PetscQuadratureDestroy(&fvm->quadrature));
1436:   fvm->quadrature = q;
1437:   PetscFunctionReturn(PETSC_SUCCESS);
1438: }

1440: /*@
1441:   PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`

1443:   Not Collective

1445:   Input Parameter:
1446: . fvm - the `PetscFV` object

1448:   Output Parameter:
1449: . q - The `PetscQuadrature`

1451:   Level: intermediate

1453: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1454: @*/
1455: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1456: {
1457:   PetscFunctionBegin;
1459:   PetscAssertPointer(q, 2);
1460:   if (!fvm->quadrature) {
1461:     /* Create default 1-point quadrature */
1462:     PetscReal *points, *weights;

1464:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature));
1465:     PetscCall(PetscCalloc1(fvm->dim, &points));
1466:     PetscCall(PetscMalloc1(1, &weights));
1467:     weights[0] = 1.0;
1468:     PetscCall(PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights));
1469:   }
1470:   *q = fvm->quadrature;
1471:   PetscFunctionReturn(PETSC_SUCCESS);
1472: }

1474: /*@
1475:   PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`

1477:   Not Collective

1479:   Input Parameters:
1480: + fvm - The `PetscFV` object
1481: - ct  - The `DMPolytopeType` for the cell

1483:   Level: intermediate

1485: .seealso: `PetscFVGetDualSpace()`, `PetscFVSetDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1486: @*/
1487: PetscErrorCode PetscFVCreateDualSpace(PetscFV fvm, DMPolytopeType ct)
1488: {
1489:   DM       K;
1490:   PetscInt dim, Nc;

1492:   PetscFunctionBegin;
1493:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1494:   PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1495:   PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1496:   PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1497:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1498:   PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1499:   PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1500:   PetscCall(DMDestroy(&K));
1501:   PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1502:   // Should we be using PetscFVGetQuadrature() here?
1503:   for (PetscInt c = 0; c < Nc; ++c) {
1504:     PetscQuadrature qc;
1505:     PetscReal      *points, *weights;

1507:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1508:     PetscCall(PetscCalloc1(dim, &points));
1509:     PetscCall(PetscCalloc1(Nc, &weights));
1510:     weights[c] = 1.0;
1511:     PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1512:     PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1513:     PetscCall(PetscQuadratureDestroy(&qc));
1514:   }
1515:   PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1516:   PetscFunctionReturn(PETSC_SUCCESS);
1517: }

1519: /*@
1520:   PetscFVGetDualSpace - Returns the `PetscDualSpace` used to define the inner product on a `PetscFV`

1522:   Not Collective

1524:   Input Parameter:
1525: . fvm - The `PetscFV` object

1527:   Output Parameter:
1528: . sp - The `PetscDualSpace` object

1530:   Level: intermediate

1532:   Developer Notes:
1533:   There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class

1535: .seealso: `PetscFVSetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1536: @*/
1537: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1538: {
1539:   PetscFunctionBegin;
1541:   PetscAssertPointer(sp, 2);
1542:   if (!fvm->dualSpace) {
1543:     PetscInt dim;

1545:     PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1546:     PetscCall(PetscFVCreateDualSpace(fvm, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE)));
1547:   }
1548:   *sp = fvm->dualSpace;
1549:   PetscFunctionReturn(PETSC_SUCCESS);
1550: }

1552: /*@
1553:   PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product

1555:   Not Collective

1557:   Input Parameters:
1558: + fvm - The `PetscFV` object
1559: - sp  - The `PetscDualSpace` object

1561:   Level: intermediate

1563:   Note:
1564:   A simple dual space is provided automatically, and the user typically will not need to override it.

1566: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1567: @*/
1568: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1569: {
1570:   PetscFunctionBegin;
1573:   PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1574:   fvm->dualSpace = sp;
1575:   PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1576:   PetscFunctionReturn(PETSC_SUCCESS);
1577: }

1579: /*@C
1580:   PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points

1582:   Not Collective

1584:   Input Parameter:
1585: . fvm - The `PetscFV` object

1587:   Output Parameter:
1588: . T - The basis function values and derivatives at quadrature points

1590:   Level: intermediate

1592:   Note:
1593: .vb
1594:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1595:   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
1596:   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
1597: .ve

1599: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1600: @*/
1601: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1602: {
1603:   PetscInt         npoints;
1604:   const PetscReal *points;

1606:   PetscFunctionBegin;
1608:   PetscAssertPointer(T, 2);
1609:   PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1610:   if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1611:   *T = fvm->T;
1612:   PetscFunctionReturn(PETSC_SUCCESS);
1613: }

1615: /*@C
1616:   PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

1618:   Not Collective

1620:   Input Parameters:
1621: + fvm     - The `PetscFV` object
1622: . nrepl   - The number of replicas
1623: . npoints - The number of tabulation points in a replica
1624: . points  - The tabulation point coordinates
1625: - K       - The order of derivative to tabulate

1627:   Output Parameter:
1628: . T - The basis function values and derivative at tabulation points

1630:   Level: intermediate

1632:   Note:
1633: .vb
1634:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1635:   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
1636:   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
1637: .ve

1639: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1640: @*/
1641: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1642: {
1643:   PetscInt pdim; // Dimension of approximation space P
1644:   PetscInt cdim; // Spatial dimension
1645:   PetscInt Nc;   // Field components
1646:   PetscInt k, p, d, c, e;

1648:   PetscFunctionBegin;
1649:   if (!npoints || K < 0) {
1650:     *T = NULL;
1651:     PetscFunctionReturn(PETSC_SUCCESS);
1652:   }
1654:   PetscAssertPointer(points, 4);
1655:   PetscAssertPointer(T, 6);
1656:   PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1657:   PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1658:   pdim = Nc;
1659:   PetscCall(PetscMalloc1(1, T));
1660:   (*T)->K    = !cdim ? 0 : K;
1661:   (*T)->Nr   = nrepl;
1662:   (*T)->Np   = npoints;
1663:   (*T)->Nb   = pdim;
1664:   (*T)->Nc   = Nc;
1665:   (*T)->cdim = cdim;
1666:   PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1667:   for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1668:   if (K >= 0) {
1669:     for (p = 0; p < nrepl * npoints; ++p)
1670:       for (d = 0; d < pdim; ++d)
1671:         for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.;
1672:   }
1673:   if (K >= 1) {
1674:     for (p = 0; p < nrepl * npoints; ++p)
1675:       for (d = 0; d < pdim; ++d)
1676:         for (c = 0; c < Nc; ++c)
1677:           for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1678:   }
1679:   if (K >= 2) {
1680:     for (p = 0; p < nrepl * npoints; ++p)
1681:       for (d = 0; d < pdim; ++d)
1682:         for (c = 0; c < Nc; ++c)
1683:           for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1684:   }
1685:   PetscFunctionReturn(PETSC_SUCCESS);
1686: }

1688: /*@C
1689:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

1691:   Input Parameters:
1692: + fvm      - The `PetscFV` object
1693: . numFaces - The number of cell faces which are not constrained
1694: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1696:   Output Parameter:
1697: . grad - the gradient

1699:   Level: advanced

1701: .seealso: `PetscFV`, `PetscFVCreate()`
1702: @*/
1703: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1704: {
1705:   PetscFunctionBegin;
1707:   PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1708:   PetscFunctionReturn(PETSC_SUCCESS);
1709: }

1711: /*@C
1712:   PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration

1714:   Not Collective

1716:   Input Parameters:
1717: + fvm         - The `PetscFV` object for the field being integrated
1718: . prob        - The `PetscDS` specifying the discretizations and continuum functions
1719: . field       - The field being integrated
1720: . Nf          - The number of faces in the chunk
1721: . fgeom       - The face geometry for each face in the chunk
1722: . neighborVol - The volume for each pair of cells in the chunk
1723: . uL          - The state from the cell on the left
1724: - uR          - The state from the cell on the right

1726:   Output Parameters:
1727: + fluxL - the left fluxes for each face
1728: - fluxR - the right fluxes for each face

1730:   Level: developer

1732: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1733: @*/
1734: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1735: {
1736:   PetscFunctionBegin;
1738:   PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1739:   PetscFunctionReturn(PETSC_SUCCESS);
1740: }

1742: /*@
1743:   PetscFVClone - Create a shallow copy of a `PetscFV` object that just references the internal objects.

1745:   Input Parameter:
1746: . fv - The initial `PetscFV`

1748:   Output Parameter:
1749: . fvNew - A clone of the `PetscFV`

1751:   Level: advanced

1753:   Notes:
1754:   This is typically used to change the number of components.

1756: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1757: @*/
1758: PetscErrorCode PetscFVClone(PetscFV fv, PetscFV *fvNew)
1759: {
1760:   PetscDualSpace  Q;
1761:   DM              K;
1762:   PetscQuadrature q;
1763:   PetscInt        Nc, cdim;

1765:   PetscFunctionBegin;
1766:   PetscCall(PetscFVGetDualSpace(fv, &Q));
1767:   PetscCall(PetscFVGetQuadrature(fv, &q));
1768:   PetscCall(PetscDualSpaceGetDM(Q, &K));

1770:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvNew));
1771:   PetscCall(PetscFVSetDualSpace(*fvNew, Q));
1772:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
1773:   PetscCall(PetscFVSetNumComponents(*fvNew, Nc));
1774:   PetscCall(PetscFVGetSpatialDimension(fv, &cdim));
1775:   PetscCall(PetscFVSetSpatialDimension(*fvNew, cdim));
1776:   PetscCall(PetscFVSetQuadrature(*fvNew, q));
1777:   PetscFunctionReturn(PETSC_SUCCESS);
1778: }

1780: /*@
1781:   PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1782:   smaller copies.

1784:   Input Parameter:
1785: . fv - The initial `PetscFV`

1787:   Output Parameter:
1788: . fvRef - The refined `PetscFV`

1790:   Level: advanced

1792:   Notes:
1793:   This is typically used to generate a preconditioner for a high order method from a lower order method on a
1794:   refined mesh having the same number of dofs (but more sparsity). It is also used to create an
1795:   interpolation between regularly refined meshes.

1797: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1798: @*/
1799: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1800: {
1801:   PetscDualSpace  Q, Qref;
1802:   DM              K, Kref;
1803:   PetscQuadrature q, qref;
1804:   DMPolytopeType  ct;
1805:   DMPlexTransform tr;
1806:   PetscReal      *v0;
1807:   PetscReal      *jac, *invjac;
1808:   PetscInt        numComp, numSubelements, s;

1810:   PetscFunctionBegin;
1811:   PetscCall(PetscFVGetDualSpace(fv, &Q));
1812:   PetscCall(PetscFVGetQuadrature(fv, &q));
1813:   PetscCall(PetscDualSpaceGetDM(Q, &K));
1814:   /* Create dual space */
1815:   PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1816:   PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1817:   PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1818:   PetscCall(DMDestroy(&Kref));
1819:   PetscCall(PetscDualSpaceSetUp(Qref));
1820:   /* Create volume */
1821:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1822:   PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1823:   PetscCall(PetscFVGetNumComponents(fv, &numComp));
1824:   PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1825:   PetscCall(PetscFVSetUp(*fvRef));
1826:   /* Create quadrature */
1827:   PetscCall(DMPlexGetCellType(K, 0, &ct));
1828:   PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1829:   PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1830:   PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1831:   PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1832:   PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1833:   for (s = 0; s < numSubelements; ++s) {
1834:     PetscQuadrature  qs;
1835:     const PetscReal *points, *weights;
1836:     PetscReal       *p, *w;
1837:     PetscInt         dim, Nc, npoints, np;

1839:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1840:     PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1841:     np = npoints / numSubelements;
1842:     PetscCall(PetscMalloc1(np * dim, &p));
1843:     PetscCall(PetscMalloc1(np * Nc, &w));
1844:     PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1845:     PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1846:     PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1847:     PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1848:     PetscCall(PetscQuadratureDestroy(&qs));
1849:   }
1850:   PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1851:   PetscCall(DMPlexTransformDestroy(&tr));
1852:   PetscCall(PetscQuadratureDestroy(&qref));
1853:   PetscCall(PetscDualSpaceDestroy(&Qref));
1854:   PetscFunctionReturn(PETSC_SUCCESS);
1855: }

1857: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1858: {
1859:   PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;

1861:   PetscFunctionBegin;
1862:   PetscCall(PetscFree(b));
1863:   PetscFunctionReturn(PETSC_SUCCESS);
1864: }

1866: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1867: {
1868:   PetscInt          Nc, c;
1869:   PetscViewerFormat format;

1871:   PetscFunctionBegin;
1872:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
1873:   PetscCall(PetscViewerGetFormat(viewer, &format));
1874:   PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1875:   PetscCall(PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc));
1876:   for (c = 0; c < Nc; c++) {
1877:     if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1878:   }
1879:   PetscFunctionReturn(PETSC_SUCCESS);
1880: }

1882: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1883: {
1884:   PetscBool iascii;

1886:   PetscFunctionBegin;
1889:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1890:   if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1891:   PetscFunctionReturn(PETSC_SUCCESS);
1892: }

1894: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1895: {
1896:   PetscFunctionBegin;
1897:   PetscFunctionReturn(PETSC_SUCCESS);
1898: }

1900: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1901: {
1902:   PetscInt dim;

1904:   PetscFunctionBegin;
1905:   PetscCall(PetscFVGetSpatialDimension(fv, &dim));
1906:   for (PetscInt f = 0; f < numFaces; ++f) {
1907:     for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
1908:   }
1909:   PetscFunctionReturn(PETSC_SUCCESS);
1910: }

1912: /*
1913:   neighborVol[f*2+0] contains the left  geom
1914:   neighborVol[f*2+1] contains the right geom
1915: */
1916: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1917: {
1918:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1919:   void              *rctx;
1920:   PetscScalar       *flux = fvm->fluxWork;
1921:   const PetscScalar *constants;
1922:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;

1924:   PetscFunctionBegin;
1925:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1926:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1927:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1928:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1929:   PetscCall(PetscDSGetContext(prob, field, &rctx));
1930:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1931:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1932:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1933:   for (f = 0; f < Nf; ++f) {
1934:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1935:     for (d = 0; d < pdim; ++d) {
1936:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1937:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1938:     }
1939:   }
1940:   PetscFunctionReturn(PETSC_SUCCESS);
1941: }

1943: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1944: {
1945:   PetscFunctionBegin;
1946:   fvm->ops->setfromoptions       = NULL;
1947:   fvm->ops->setup                = PetscFVSetUp_Upwind;
1948:   fvm->ops->view                 = PetscFVView_Upwind;
1949:   fvm->ops->destroy              = PetscFVDestroy_Upwind;
1950:   fvm->ops->computegradient      = PetscFVComputeGradient_Upwind;
1951:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1952:   PetscFunctionReturn(PETSC_SUCCESS);
1953: }

1955: /*MC
1956:   PETSCFVUPWIND = "upwind" - A `PetscFV` implementation

1958:   Level: intermediate

1960: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1961: M*/

1963: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1964: {
1965:   PetscFV_Upwind *b;

1967:   PetscFunctionBegin;
1969:   PetscCall(PetscNew(&b));
1970:   fvm->data = b;

1972:   PetscCall(PetscFVInitialize_Upwind(fvm));
1973:   PetscFunctionReturn(PETSC_SUCCESS);
1974: }

1976: #include <petscblaslapack.h>

1978: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1979: {
1980:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;

1982:   PetscFunctionBegin;
1983:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1984:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1985:   PetscCall(PetscFree(ls));
1986:   PetscFunctionReturn(PETSC_SUCCESS);
1987: }

1989: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1990: {
1991:   PetscInt          Nc, c;
1992:   PetscViewerFormat format;

1994:   PetscFunctionBegin;
1995:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
1996:   PetscCall(PetscViewerGetFormat(viewer, &format));
1997:   PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1998:   PetscCall(PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc));
1999:   for (c = 0; c < Nc; c++) {
2000:     if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
2001:   }
2002:   PetscFunctionReturn(PETSC_SUCCESS);
2003: }

2005: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2006: {
2007:   PetscBool iascii;

2009:   PetscFunctionBegin;
2012:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2013:   if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
2014:   PetscFunctionReturn(PETSC_SUCCESS);
2015: }

2017: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2018: {
2019:   PetscFunctionBegin;
2020:   PetscFunctionReturn(PETSC_SUCCESS);
2021: }

2023: /* Overwrites A. Can only handle full-rank problems with m>=n */
2024: static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2025: {
2026:   PetscBool    debug = PETSC_FALSE;
2027:   PetscBLASInt M, N, K, lda, ldb, ldwork, info;
2028:   PetscScalar *R, *Q, *Aback, Alpha;

2030:   PetscFunctionBegin;
2031:   if (debug) {
2032:     PetscCall(PetscMalloc1(m * n, &Aback));
2033:     PetscCall(PetscArraycpy(Aback, A, m * n));
2034:   }

2036:   PetscCall(PetscBLASIntCast(m, &M));
2037:   PetscCall(PetscBLASIntCast(n, &N));
2038:   PetscCall(PetscBLASIntCast(mstride, &lda));
2039:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
2040:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2041:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
2042:   PetscCall(PetscFPTrapPop());
2043:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
2044:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2046:   /* Extract an explicit representation of Q */
2047:   Q = Ainv;
2048:   PetscCall(PetscArraycpy(Q, A, mstride * n));
2049:   K = N; /* full rank */
2050:   PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
2051:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");

2053:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2054:   Alpha = 1.0;
2055:   ldb   = lda;
2056:   BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
2057:   /* Ainv is Q, overwritten with inverse */

2059:   if (debug) { /* Check that pseudo-inverse worked */
2060:     PetscScalar  Beta = 0.0;
2061:     PetscBLASInt ldc;
2062:     K   = N;
2063:     ldc = N;
2064:     BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
2065:     PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
2066:     PetscCall(PetscFree(Aback));
2067:   }
2068:   PetscFunctionReturn(PETSC_SUCCESS);
2069: }

2071: /* Overwrites A. Can handle degenerate problems and m<n. */
2072: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
2073: {
2074:   PetscScalar *Brhs;
2075:   PetscScalar *tmpwork;
2076:   PetscReal    rcond;
2077: #if defined(PETSC_USE_COMPLEX)
2078:   PetscInt   rworkSize;
2079:   PetscReal *rwork;
2080: #endif
2081:   PetscInt     i, j, maxmn;
2082:   PetscBLASInt M, N, lda, ldb, ldwork;
2083:   PetscBLASInt nrhs, irank, info;

2085:   PetscFunctionBegin;
2086:   /* initialize to identity */
2087:   tmpwork = work;
2088:   Brhs    = Ainv;
2089:   maxmn   = PetscMax(m, n);
2090:   for (j = 0; j < maxmn; j++) {
2091:     for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2092:   }

2094:   PetscCall(PetscBLASIntCast(m, &M));
2095:   PetscCall(PetscBLASIntCast(n, &N));
2096:   PetscCall(PetscBLASIntCast(mstride, &lda));
2097:   PetscCall(PetscBLASIntCast(maxmn, &ldb));
2098:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
2099:   rcond = -1;
2100:   nrhs  = M;
2101: #if defined(PETSC_USE_COMPLEX)
2102:   rworkSize = 5 * PetscMin(M, N);
2103:   PetscCall(PetscMalloc1(rworkSize, &rwork));
2104:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2105:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2106:   PetscCall(PetscFPTrapPop());
2107:   PetscCall(PetscFree(rwork));
2108: #else
2109:   nrhs = M;
2110:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2111:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2112:   PetscCall(PetscFPTrapPop());
2113: #endif
2114:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2115:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2116:   PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2117:   PetscFunctionReturn(PETSC_SUCCESS);
2118: }

2120: #if 0
2121: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2122: {
2123:   PetscReal       grad[2] = {0, 0};
2124:   const PetscInt *faces;
2125:   PetscInt        numFaces, f;

2127:   PetscFunctionBegin;
2128:   PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2129:   PetscCall(DMPlexGetCone(dm, cell, &faces));
2130:   for (f = 0; f < numFaces; ++f) {
2131:     const PetscInt *fcells;
2132:     const CellGeom *cg1;
2133:     const FaceGeom *fg;

2135:     PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2136:     PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2137:     for (i = 0; i < 2; ++i) {
2138:       PetscScalar du;

2140:       if (fcells[i] == c) continue;
2141:       PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2142:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2143:       grad[0] += fg->grad[!i][0] * du;
2144:       grad[1] += fg->grad[!i][1] * du;
2145:     }
2146:   }
2147:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2148:   PetscFunctionReturn(PETSC_SUCCESS);
2149: }
2150: #endif

2152: /*
2153:   PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell

2155:   Input Parameters:
2156: + fvm      - The `PetscFV` object
2157: . numFaces - The number of cell faces which are not constrained
2158: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

2160:   Level: developer

2162: .seealso: `PetscFV`, `PetscFVCreate()`
2163: */
2164: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2165: {
2166:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *)fvm->data;
2167:   const PetscBool       useSVD   = PETSC_TRUE;
2168:   const PetscInt        maxFaces = ls->maxFaces;
2169:   PetscInt              dim, f, d;

2171:   PetscFunctionBegin;
2172:   if (numFaces > maxFaces) {
2173:     PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2174:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2175:   }
2176:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2177:   for (f = 0; f < numFaces; ++f) {
2178:     for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2179:   }
2180:   /* Overwrites B with garbage, returns Binv in row-major format */
2181:   if (useSVD) {
2182:     PetscInt maxmn = PetscMax(numFaces, dim);
2183:     PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2184:     /* Binv shaped in column-major, coldim=maxmn.*/
2185:     for (f = 0; f < numFaces; ++f) {
2186:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2187:     }
2188:   } else {
2189:     PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2190:     /* Binv shaped in row-major, rowdim=maxFaces.*/
2191:     for (f = 0; f < numFaces; ++f) {
2192:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2193:     }
2194:   }
2195:   PetscFunctionReturn(PETSC_SUCCESS);
2196: }

2198: /*
2199:   neighborVol[f*2+0] contains the left  geom
2200:   neighborVol[f*2+1] contains the right geom
2201: */
2202: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2203: {
2204:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2205:   void              *rctx;
2206:   PetscScalar       *flux = fvm->fluxWork;
2207:   const PetscScalar *constants;
2208:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;

2210:   PetscFunctionBegin;
2211:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2212:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2213:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2214:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2215:   PetscCall(PetscDSGetContext(prob, field, &rctx));
2216:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2217:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2218:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2219:   for (f = 0; f < Nf; ++f) {
2220:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2221:     for (d = 0; d < pdim; ++d) {
2222:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2223:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2224:     }
2225:   }
2226:   PetscFunctionReturn(PETSC_SUCCESS);
2227: }

2229: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2230: {
2231:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2232:   PetscInt              dim, m, n, nrhs, minmn, maxmn;

2234:   PetscFunctionBegin;
2236:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2237:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2238:   ls->maxFaces = maxFaces;
2239:   m            = ls->maxFaces;
2240:   n            = dim;
2241:   nrhs         = ls->maxFaces;
2242:   minmn        = PetscMin(m, n);
2243:   maxmn        = PetscMax(m, n);
2244:   ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2245:   PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2246:   PetscFunctionReturn(PETSC_SUCCESS);
2247: }

2249: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2250: {
2251:   PetscFunctionBegin;
2252:   fvm->ops->setfromoptions       = NULL;
2253:   fvm->ops->setup                = PetscFVSetUp_LeastSquares;
2254:   fvm->ops->view                 = PetscFVView_LeastSquares;
2255:   fvm->ops->destroy              = PetscFVDestroy_LeastSquares;
2256:   fvm->ops->computegradient      = PetscFVComputeGradient_LeastSquares;
2257:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2258:   PetscFunctionReturn(PETSC_SUCCESS);
2259: }

2261: /*MC
2262:   PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation

2264:   Level: intermediate

2266: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2267: M*/

2269: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2270: {
2271:   PetscFV_LeastSquares *ls;

2273:   PetscFunctionBegin;
2275:   PetscCall(PetscNew(&ls));
2276:   fvm->data = ls;

2278:   ls->maxFaces = -1;
2279:   ls->workSize = -1;
2280:   ls->B        = NULL;
2281:   ls->Binv     = NULL;
2282:   ls->tau      = NULL;
2283:   ls->work     = NULL;

2285:   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2286:   PetscCall(PetscFVInitialize_LeastSquares(fvm));
2287:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2288:   PetscFunctionReturn(PETSC_SUCCESS);
2289: }

2291: /*@
2292:   PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction

2294:   Not Collective

2296:   Input Parameters:
2297: + fvm      - The `PetscFV` object
2298: - maxFaces - The maximum number of cell faces

2300:   Level: intermediate

2302: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2303: @*/
2304: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2305: {
2306:   PetscFunctionBegin;
2308:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2309:   PetscFunctionReturn(PETSC_SUCCESS);
2310: }