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:   if (lim->ops->destroy) {
 88:     PetscUseTypeMethod(lim, destroy);
 89:     lim->ops->destroy = NULL;
 90:   }
 91:   PetscCall((*r)(lim));
 92:   PetscCall(PetscObjectChangeTypeName((PetscObject)lim, name));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

 99:   Not Collective

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

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

107:   Level: intermediate

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

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

124:   Collective

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

131:   Level: intermediate

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

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

146:   Collective

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

152:   Level: beginner

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

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

168:   Collective

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

173:   Level: intermediate

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

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

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

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

207:   Collective

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

212:   Level: intermediate

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

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

227:   Collective

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

232:   Level: beginner

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

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

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

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

256:   Collective

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

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

264:   Level: beginner

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

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

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

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

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

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

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

294:   Level: beginner

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

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

303:  The second order TVD region is bounded by

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

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

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

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

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

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

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

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

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

324:  becomes

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

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

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

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

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

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

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

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

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

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

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

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

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

394:   Level: intermediate

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

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

403:   PetscFunctionBegin;
405:   PetscCall(PetscNew(&l));
406:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

462:   Level: intermediate

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

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

471:   PetscFunctionBegin;
473:   PetscCall(PetscNew(&l));
474:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

530:   Level: intermediate

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

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

539:   PetscFunctionBegin;
541:   PetscCall(PetscNew(&l));
542:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

598:   Level: intermediate

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

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

607:   PetscFunctionBegin;
609:   PetscCall(PetscNew(&l));
610:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

666:   Level: intermediate

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

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

675:   PetscFunctionBegin;
677:   PetscCall(PetscNew(&l));
678:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

734:   Level: intermediate

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

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

743:   PetscFunctionBegin;
745:   PetscCall(PetscNew(&l));
746:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

802:   Level: intermediate

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

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

811:   PetscFunctionBegin;
813:   PetscCall(PetscNew(&l));
814:   lim->data = l;

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

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

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

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

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

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

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

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

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

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

871:   Level: intermediate

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

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

880:   PetscFunctionBegin;
882:   PetscCall(PetscNew(&l));
883:   lim->data = l;

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

889: PetscClassId PETSCFV_CLASSID = 0;

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

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

897:   Not Collective

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

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

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

918:   Level: advanced

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

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

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

935:   Collective

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

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

944:   Level: intermediate

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

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

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

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

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

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

973:   Not Collective

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

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

981:   Level: intermediate

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

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

998:   Collective

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

1005:   Level: intermediate

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

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

1020:   Collective

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

1026:   Level: beginner

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

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

1042:   Collective

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

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

1050:   Level: intermediate

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

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

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

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

1086:   Collective

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

1091:   Level: intermediate

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

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

1107:   Collective

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

1112:   Level: beginner

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

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

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

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

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

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

1146:   Collective

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

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

1154:   Level: beginner

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

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

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

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

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

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

1184:   Logically Collective

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

1190:   Level: intermediate

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

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

1208:   Not Collective

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

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

1216:   Level: intermediate

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

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

1232:   Logically Collective

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

1238:   Level: intermediate

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

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

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

1262:   Not Collective

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

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

1270:   Level: intermediate

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

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

1286:   Logically Collective

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

1293:   Level: intermediate

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

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

1308:   Logically Collective

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

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

1317:   Level: intermediate

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

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

1331:   Logically Collective

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

1337:   Level: intermediate

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

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

1352:   Not Collective

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

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

1360:   Level: intermediate

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

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

1376:   Logically Collective

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

1382:   Level: intermediate

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

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

1397:   Not Collective

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

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

1405:   Level: intermediate

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

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

1421:   Logically Collective

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

1427:   Level: intermediate

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

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

1444:   Not Collective

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

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

1452:   Level: intermediate

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

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

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

1478:   Not Collective

1480:   Input Parameter:
1481: . fvm - The `PetscFV` object

1483:   Output Parameter:
1484: . sp - The `PetscDualSpace` object

1486:   Level: intermediate

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

1491: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1492: @*/
1493: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1494: {
1495:   PetscFunctionBegin;
1497:   PetscAssertPointer(sp, 2);
1498:   if (!fvm->dualSpace) {
1499:     DM       K;
1500:     PetscInt dim, Nc, c;

1502:     PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1503:     PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1504:     PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace));
1505:     PetscCall(PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE));
1506:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K));
1507:     PetscCall(PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc));
1508:     PetscCall(PetscDualSpaceSetDM(fvm->dualSpace, K));
1509:     PetscCall(DMDestroy(&K));
1510:     PetscCall(PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc));
1511:     /* Should we be using PetscFVGetQuadrature() here? */
1512:     for (c = 0; c < Nc; ++c) {
1513:       PetscQuadrature qc;
1514:       PetscReal      *points, *weights;

1516:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qc));
1517:       PetscCall(PetscCalloc1(dim, &points));
1518:       PetscCall(PetscCalloc1(Nc, &weights));
1519:       weights[c] = 1.0;
1520:       PetscCall(PetscQuadratureSetData(qc, dim, Nc, 1, points, weights));
1521:       PetscCall(PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc));
1522:       PetscCall(PetscQuadratureDestroy(&qc));
1523:     }
1524:     PetscCall(PetscDualSpaceSetUp(fvm->dualSpace));
1525:   }
1526:   *sp = fvm->dualSpace;
1527:   PetscFunctionReturn(PETSC_SUCCESS);
1528: }

1530: /*@
1531:   PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product

1533:   Not Collective

1535:   Input Parameters:
1536: + fvm - The `PetscFV` object
1537: - sp  - The `PetscDualSpace` object

1539:   Level: intermediate

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

1544: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1545: @*/
1546: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1547: {
1548:   PetscFunctionBegin;
1551:   PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
1552:   fvm->dualSpace = sp;
1553:   PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
1554:   PetscFunctionReturn(PETSC_SUCCESS);
1555: }

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

1560:   Not Collective

1562:   Input Parameter:
1563: . fvm - The `PetscFV` object

1565:   Output Parameter:
1566: . T - The basis function values and derivatives at quadrature points

1568:   Level: intermediate

1570:   Note:
1571: .vb
1572:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1573:   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
1574:   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
1575: .ve

1577: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1578: @*/
1579: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1580: {
1581:   PetscInt         npoints;
1582:   const PetscReal *points;

1584:   PetscFunctionBegin;
1586:   PetscAssertPointer(T, 2);
1587:   PetscCall(PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL));
1588:   if (!fvm->T) PetscCall(PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T));
1589:   *T = fvm->T;
1590:   PetscFunctionReturn(PETSC_SUCCESS);
1591: }

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

1596:   Not Collective

1598:   Input Parameters:
1599: + fvm     - The `PetscFV` object
1600: . nrepl   - The number of replicas
1601: . npoints - The number of tabulation points in a replica
1602: . points  - The tabulation point coordinates
1603: - K       - The order of derivative to tabulate

1605:   Output Parameter:
1606: . T - The basis function values and derivative at tabulation points

1608:   Level: intermediate

1610:   Note:
1611: .vb
1612:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1613:   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
1614:   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
1615: .ve

1617: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1618: @*/
1619: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1620: {
1621:   PetscInt pdim = 1; /* Dimension of approximation space P */
1622:   PetscInt cdim;     /* Spatial dimension */
1623:   PetscInt Nc;       /* Field components */
1624:   PetscInt k, p, d, c, e;

1626:   PetscFunctionBegin;
1627:   if (!npoints || K < 0) {
1628:     *T = NULL;
1629:     PetscFunctionReturn(PETSC_SUCCESS);
1630:   }
1632:   PetscAssertPointer(points, 4);
1633:   PetscAssertPointer(T, 6);
1634:   PetscCall(PetscFVGetSpatialDimension(fvm, &cdim));
1635:   PetscCall(PetscFVGetNumComponents(fvm, &Nc));
1636:   PetscCall(PetscMalloc1(1, T));
1637:   (*T)->K    = !cdim ? 0 : K;
1638:   (*T)->Nr   = nrepl;
1639:   (*T)->Np   = npoints;
1640:   (*T)->Nb   = pdim;
1641:   (*T)->Nc   = Nc;
1642:   (*T)->cdim = cdim;
1643:   PetscCall(PetscMalloc1((*T)->K + 1, &(*T)->T));
1644:   for (k = 0; k <= (*T)->K; ++k) PetscCall(PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]));
1645:   if (K >= 0) {
1646:     for (p = 0; p < nrepl * npoints; ++p)
1647:       for (d = 0; d < pdim; ++d)
1648:         for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1649:   }
1650:   if (K >= 1) {
1651:     for (p = 0; p < nrepl * npoints; ++p)
1652:       for (d = 0; d < pdim; ++d)
1653:         for (c = 0; c < Nc; ++c)
1654:           for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1655:   }
1656:   if (K >= 2) {
1657:     for (p = 0; p < nrepl * npoints; ++p)
1658:       for (d = 0; d < pdim; ++d)
1659:         for (c = 0; c < Nc; ++c)
1660:           for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1661:   }
1662:   PetscFunctionReturn(PETSC_SUCCESS);
1663: }

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

1668:   Input Parameters:
1669: + fvm      - The `PetscFV` object
1670: . numFaces - The number of cell faces which are not constrained
1671: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1673:   Output Parameter:
1674: . grad - the gradient

1676:   Level: advanced

1678: .seealso: `PetscFV`, `PetscFVCreate()`
1679: @*/
1680: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1681: {
1682:   PetscFunctionBegin;
1684:   PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1685:   PetscFunctionReturn(PETSC_SUCCESS);
1686: }

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

1691:   Not Collective

1693:   Input Parameters:
1694: + fvm         - The `PetscFV` object for the field being integrated
1695: . prob        - The `PetscDS` specifying the discretizations and continuum functions
1696: . field       - The field being integrated
1697: . Nf          - The number of faces in the chunk
1698: . fgeom       - The face geometry for each face in the chunk
1699: . neighborVol - The volume for each pair of cells in the chunk
1700: . uL          - The state from the cell on the left
1701: - uR          - The state from the cell on the right

1703:   Output Parameters:
1704: + fluxL - the left fluxes for each face
1705: - fluxR - the right fluxes for each face

1707:   Level: developer

1709: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1710: @*/
1711: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1712: {
1713:   PetscFunctionBegin;
1715:   PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1716:   PetscFunctionReturn(PETSC_SUCCESS);
1717: }

1719: /*@
1720:   PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
1721:   smaller copies.

1723:   Input Parameter:
1724: . fv - The initial `PetscFV`

1726:   Output Parameter:
1727: . fvRef - The refined `PetscFV`

1729:   Level: advanced

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

1736: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1737: @*/
1738: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1739: {
1740:   PetscDualSpace  Q, Qref;
1741:   DM              K, Kref;
1742:   PetscQuadrature q, qref;
1743:   DMPolytopeType  ct;
1744:   DMPlexTransform tr;
1745:   PetscReal      *v0;
1746:   PetscReal      *jac, *invjac;
1747:   PetscInt        numComp, numSubelements, s;

1749:   PetscFunctionBegin;
1750:   PetscCall(PetscFVGetDualSpace(fv, &Q));
1751:   PetscCall(PetscFVGetQuadrature(fv, &q));
1752:   PetscCall(PetscDualSpaceGetDM(Q, &K));
1753:   /* Create dual space */
1754:   PetscCall(PetscDualSpaceDuplicate(Q, &Qref));
1755:   PetscCall(DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref));
1756:   PetscCall(PetscDualSpaceSetDM(Qref, Kref));
1757:   PetscCall(DMDestroy(&Kref));
1758:   PetscCall(PetscDualSpaceSetUp(Qref));
1759:   /* Create volume */
1760:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef));
1761:   PetscCall(PetscFVSetDualSpace(*fvRef, Qref));
1762:   PetscCall(PetscFVGetNumComponents(fv, &numComp));
1763:   PetscCall(PetscFVSetNumComponents(*fvRef, numComp));
1764:   PetscCall(PetscFVSetUp(*fvRef));
1765:   /* Create quadrature */
1766:   PetscCall(DMPlexGetCellType(K, 0, &ct));
1767:   PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
1768:   PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
1769:   PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac));
1770:   PetscCall(PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref));
1771:   PetscCall(PetscDualSpaceSimpleSetDimension(Qref, numSubelements));
1772:   for (s = 0; s < numSubelements; ++s) {
1773:     PetscQuadrature  qs;
1774:     const PetscReal *points, *weights;
1775:     PetscReal       *p, *w;
1776:     PetscInt         dim, Nc, npoints, np;

1778:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qs));
1779:     PetscCall(PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights));
1780:     np = npoints / numSubelements;
1781:     PetscCall(PetscMalloc1(np * dim, &p));
1782:     PetscCall(PetscMalloc1(np * Nc, &w));
1783:     PetscCall(PetscArraycpy(p, &points[s * np * dim], np * dim));
1784:     PetscCall(PetscArraycpy(w, &weights[s * np * Nc], np * Nc));
1785:     PetscCall(PetscQuadratureSetData(qs, dim, Nc, np, p, w));
1786:     PetscCall(PetscDualSpaceSimpleSetFunctional(Qref, s, qs));
1787:     PetscCall(PetscQuadratureDestroy(&qs));
1788:   }
1789:   PetscCall(PetscFVSetQuadrature(*fvRef, qref));
1790:   PetscCall(DMPlexTransformDestroy(&tr));
1791:   PetscCall(PetscQuadratureDestroy(&qref));
1792:   PetscCall(PetscDualSpaceDestroy(&Qref));
1793:   PetscFunctionReturn(PETSC_SUCCESS);
1794: }

1796: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1797: {
1798:   PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;

1800:   PetscFunctionBegin;
1801:   PetscCall(PetscFree(b));
1802:   PetscFunctionReturn(PETSC_SUCCESS);
1803: }

1805: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1806: {
1807:   PetscInt          Nc, c;
1808:   PetscViewerFormat format;

1810:   PetscFunctionBegin;
1811:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
1812:   PetscCall(PetscViewerGetFormat(viewer, &format));
1813:   PetscCall(PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n"));
1814:   PetscCall(PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc));
1815:   for (c = 0; c < Nc; c++) {
1816:     if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1817:   }
1818:   PetscFunctionReturn(PETSC_SUCCESS);
1819: }

1821: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1822: {
1823:   PetscBool iascii;

1825:   PetscFunctionBegin;
1828:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1829:   if (iascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
1830:   PetscFunctionReturn(PETSC_SUCCESS);
1831: }

1833: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1834: {
1835:   PetscFunctionBegin;
1836:   PetscFunctionReturn(PETSC_SUCCESS);
1837: }

1839: /*
1840:   neighborVol[f*2+0] contains the left  geom
1841:   neighborVol[f*2+1] contains the right geom
1842: */
1843: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1844: {
1845:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1846:   void              *rctx;
1847:   PetscScalar       *flux = fvm->fluxWork;
1848:   const PetscScalar *constants;
1849:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;

1851:   PetscFunctionBegin;
1852:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
1853:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
1854:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
1855:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
1856:   PetscCall(PetscDSGetContext(prob, field, &rctx));
1857:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
1858:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
1859:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
1860:   for (f = 0; f < Nf; ++f) {
1861:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1862:     for (d = 0; d < pdim; ++d) {
1863:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1864:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1865:     }
1866:   }
1867:   PetscFunctionReturn(PETSC_SUCCESS);
1868: }

1870: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1871: {
1872:   PetscFunctionBegin;
1873:   fvm->ops->setfromoptions       = NULL;
1874:   fvm->ops->setup                = PetscFVSetUp_Upwind;
1875:   fvm->ops->view                 = PetscFVView_Upwind;
1876:   fvm->ops->destroy              = PetscFVDestroy_Upwind;
1877:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1878:   PetscFunctionReturn(PETSC_SUCCESS);
1879: }

1881: /*MC
1882:   PETSCFVUPWIND = "upwind" - A `PetscFV` implementation

1884:   Level: intermediate

1886: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1887: M*/

1889: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1890: {
1891:   PetscFV_Upwind *b;

1893:   PetscFunctionBegin;
1895:   PetscCall(PetscNew(&b));
1896:   fvm->data = b;

1898:   PetscCall(PetscFVInitialize_Upwind(fvm));
1899:   PetscFunctionReturn(PETSC_SUCCESS);
1900: }

1902: #include <petscblaslapack.h>

1904: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1905: {
1906:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;

1908:   PetscFunctionBegin;
1909:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
1910:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
1911:   PetscCall(PetscFree(ls));
1912:   PetscFunctionReturn(PETSC_SUCCESS);
1913: }

1915: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1916: {
1917:   PetscInt          Nc, c;
1918:   PetscViewerFormat format;

1920:   PetscFunctionBegin;
1921:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
1922:   PetscCall(PetscViewerGetFormat(viewer, &format));
1923:   PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
1924:   PetscCall(PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc));
1925:   for (c = 0; c < Nc; c++) {
1926:     if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
1927:   }
1928:   PetscFunctionReturn(PETSC_SUCCESS);
1929: }

1931: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1932: {
1933:   PetscBool iascii;

1935:   PetscFunctionBegin;
1938:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1939:   if (iascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
1940:   PetscFunctionReturn(PETSC_SUCCESS);
1941: }

1943: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1944: {
1945:   PetscFunctionBegin;
1946:   PetscFunctionReturn(PETSC_SUCCESS);
1947: }

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

1956:   PetscFunctionBegin;
1957:   if (debug) {
1958:     PetscCall(PetscMalloc1(m * n, &Aback));
1959:     PetscCall(PetscArraycpy(Aback, A, m * n));
1960:   }

1962:   PetscCall(PetscBLASIntCast(m, &M));
1963:   PetscCall(PetscBLASIntCast(n, &N));
1964:   PetscCall(PetscBLASIntCast(mstride, &lda));
1965:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
1966:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1967:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1968:   PetscCall(PetscFPTrapPop());
1969:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
1970:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1972:   /* Extract an explicit representation of Q */
1973:   Q = Ainv;
1974:   PetscCall(PetscArraycpy(Q, A, mstride * n));
1975:   K = N; /* full rank */
1976:   PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
1977:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");

1979:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1980:   Alpha = 1.0;
1981:   ldb   = lda;
1982:   BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1983:   /* Ainv is Q, overwritten with inverse */

1985:   if (debug) { /* Check that pseudo-inverse worked */
1986:     PetscScalar  Beta = 0.0;
1987:     PetscBLASInt ldc;
1988:     K   = N;
1989:     ldc = N;
1990:     BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1991:     PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
1992:     PetscCall(PetscFree(Aback));
1993:   }
1994:   PetscFunctionReturn(PETSC_SUCCESS);
1995: }

1997: /* Overwrites A. Can handle degenerate problems and m<n. */
1998: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1999: {
2000:   PetscScalar *Brhs;
2001:   PetscScalar *tmpwork;
2002:   PetscReal    rcond;
2003: #if defined(PETSC_USE_COMPLEX)
2004:   PetscInt   rworkSize;
2005:   PetscReal *rwork;
2006: #endif
2007:   PetscInt     i, j, maxmn;
2008:   PetscBLASInt M, N, lda, ldb, ldwork;
2009:   PetscBLASInt nrhs, irank, info;

2011:   PetscFunctionBegin;
2012:   /* initialize to identity */
2013:   tmpwork = work;
2014:   Brhs    = Ainv;
2015:   maxmn   = PetscMax(m, n);
2016:   for (j = 0; j < maxmn; j++) {
2017:     for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
2018:   }

2020:   PetscCall(PetscBLASIntCast(m, &M));
2021:   PetscCall(PetscBLASIntCast(n, &N));
2022:   PetscCall(PetscBLASIntCast(mstride, &lda));
2023:   PetscCall(PetscBLASIntCast(maxmn, &ldb));
2024:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
2025:   rcond = -1;
2026:   nrhs  = M;
2027: #if defined(PETSC_USE_COMPLEX)
2028:   rworkSize = 5 * PetscMin(M, N);
2029:   PetscCall(PetscMalloc1(rworkSize, &rwork));
2030:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2031:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
2032:   PetscCall(PetscFPTrapPop());
2033:   PetscCall(PetscFree(rwork));
2034: #else
2035:   nrhs = M;
2036:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2037:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
2038:   PetscCall(PetscFPTrapPop());
2039: #endif
2040:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
2041:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2042:   PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two colinear points");
2043:   PetscFunctionReturn(PETSC_SUCCESS);
2044: }

2046: #if 0
2047: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2048: {
2049:   PetscReal       grad[2] = {0, 0};
2050:   const PetscInt *faces;
2051:   PetscInt        numFaces, f;

2053:   PetscFunctionBegin;
2054:   PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
2055:   PetscCall(DMPlexGetCone(dm, cell, &faces));
2056:   for (f = 0; f < numFaces; ++f) {
2057:     const PetscInt *fcells;
2058:     const CellGeom *cg1;
2059:     const FaceGeom *fg;

2061:     PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2062:     PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
2063:     for (i = 0; i < 2; ++i) {
2064:       PetscScalar du;

2066:       if (fcells[i] == c) continue;
2067:       PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
2068:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2069:       grad[0] += fg->grad[!i][0] * du;
2070:       grad[1] += fg->grad[!i][1] * du;
2071:     }
2072:   }
2073:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
2074:   PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2076: #endif

2078: /*
2079:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

2081:   Input Parameters:
2082: + fvm      - The `PetscFV` object
2083: . numFaces - The number of cell faces which are not constrained
2084: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

2086:   Level: developer

2088: .seealso: `PetscFV`, `PetscFVCreate()`
2089: */
2090: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2091: {
2092:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *)fvm->data;
2093:   const PetscBool       useSVD   = PETSC_TRUE;
2094:   const PetscInt        maxFaces = ls->maxFaces;
2095:   PetscInt              dim, f, d;

2097:   PetscFunctionBegin;
2098:   if (numFaces > maxFaces) {
2099:     PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2100:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
2101:   }
2102:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2103:   for (f = 0; f < numFaces; ++f) {
2104:     for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
2105:   }
2106:   /* Overwrites B with garbage, returns Binv in row-major format */
2107:   if (useSVD) {
2108:     PetscInt maxmn = PetscMax(numFaces, dim);
2109:     PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2110:     /* Binv shaped in column-major, coldim=maxmn.*/
2111:     for (f = 0; f < numFaces; ++f) {
2112:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2113:     }
2114:   } else {
2115:     PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
2116:     /* Binv shaped in row-major, rowdim=maxFaces.*/
2117:     for (f = 0; f < numFaces; ++f) {
2118:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2119:     }
2120:   }
2121:   PetscFunctionReturn(PETSC_SUCCESS);
2122: }

2124: /*
2125:   neighborVol[f*2+0] contains the left  geom
2126:   neighborVol[f*2+1] contains the right geom
2127: */
2128: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2129: {
2130:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2131:   void              *rctx;
2132:   PetscScalar       *flux = fvm->fluxWork;
2133:   const PetscScalar *constants;
2134:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;

2136:   PetscFunctionBegin;
2137:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
2138:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
2139:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
2140:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
2141:   PetscCall(PetscDSGetContext(prob, field, &rctx));
2142:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
2143:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2144:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2145:   for (f = 0; f < Nf; ++f) {
2146:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2147:     for (d = 0; d < pdim; ++d) {
2148:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2149:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2150:     }
2151:   }
2152:   PetscFunctionReturn(PETSC_SUCCESS);
2153: }

2155: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2156: {
2157:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2158:   PetscInt              dim, m, n, nrhs, minmn, maxmn;

2160:   PetscFunctionBegin;
2162:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
2163:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
2164:   ls->maxFaces = maxFaces;
2165:   m            = ls->maxFaces;
2166:   n            = dim;
2167:   nrhs         = ls->maxFaces;
2168:   minmn        = PetscMin(m, n);
2169:   maxmn        = PetscMax(m, n);
2170:   ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2171:   PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
2172:   PetscFunctionReturn(PETSC_SUCCESS);
2173: }

2175: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2176: {
2177:   PetscFunctionBegin;
2178:   fvm->ops->setfromoptions       = NULL;
2179:   fvm->ops->setup                = PetscFVSetUp_LeastSquares;
2180:   fvm->ops->view                 = PetscFVView_LeastSquares;
2181:   fvm->ops->destroy              = PetscFVDestroy_LeastSquares;
2182:   fvm->ops->computegradient      = PetscFVComputeGradient_LeastSquares;
2183:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2184:   PetscFunctionReturn(PETSC_SUCCESS);
2185: }

2187: /*MC
2188:   PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation

2190:   Level: intermediate

2192: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2193: M*/

2195: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2196: {
2197:   PetscFV_LeastSquares *ls;

2199:   PetscFunctionBegin;
2201:   PetscCall(PetscNew(&ls));
2202:   fvm->data = ls;

2204:   ls->maxFaces = -1;
2205:   ls->workSize = -1;
2206:   ls->B        = NULL;
2207:   ls->Binv     = NULL;
2208:   ls->tau      = NULL;
2209:   ls->work     = NULL;

2211:   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
2212:   PetscCall(PetscFVInitialize_LeastSquares(fvm));
2213:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
2214:   PetscFunctionReturn(PETSC_SUCCESS);
2215: }

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

2220:   Not Collective

2222:   Input Parameters:
2223: + fvm      - The `PetscFV` object
2224: - maxFaces - The maximum number of cell faces

2226:   Level: intermediate

2228: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2229: @*/
2230: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2231: {
2232:   PetscFunctionBegin;
2234:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2235:   PetscFunctionReturn(PETSC_SUCCESS);
2236: }