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, No Fortran Support

     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: /*@
     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: /*@
     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: /*@
    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: /*@
    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: /*@
    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:   PetscCall(PetscFVInitializePackage());

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

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

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

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

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

    292:   Level: beginner

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

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

    301:  The second order TVD region is bounded by

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

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

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

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

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

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

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

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

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

    322:  becomes

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

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

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

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

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

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

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

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

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

    365:   PetscFunctionBegin;
    368:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    369:   if (isascii) PetscCall(PetscLimiterView_Sin_Ascii(lim, viewer));
    370:   PetscFunctionReturn(PETSC_SUCCESS);
    371: }

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

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

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

    392:   Level: intermediate

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

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

    401:   PetscFunctionBegin;
    403:   PetscCall(PetscNew(&l));
    404:   lim->data = l;

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

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

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

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

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

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

    433:   PetscFunctionBegin;
    436:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    437:   if (isascii) PetscCall(PetscLimiterView_Zero_Ascii(lim, viewer));
    438:   PetscFunctionReturn(PETSC_SUCCESS);
    439: }

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

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

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

    460:   Level: intermediate

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

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

    469:   PetscFunctionBegin;
    471:   PetscCall(PetscNew(&l));
    472:   lim->data = l;

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

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

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

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

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

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

    501:   PetscFunctionBegin;
    504:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    505:   if (isascii) PetscCall(PetscLimiterView_None_Ascii(lim, viewer));
    506:   PetscFunctionReturn(PETSC_SUCCESS);
    507: }

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

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

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

    528:   Level: intermediate

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

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

    537:   PetscFunctionBegin;
    539:   PetscCall(PetscNew(&l));
    540:   lim->data = l;

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

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

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

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

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

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

    569:   PetscFunctionBegin;
    572:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    573:   if (isascii) PetscCall(PetscLimiterView_Minmod_Ascii(lim, viewer));
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

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

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

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

    596:   Level: intermediate

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

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

    605:   PetscFunctionBegin;
    607:   PetscCall(PetscNew(&l));
    608:   lim->data = l;

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

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

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

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

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

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

    637:   PetscFunctionBegin;
    640:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    641:   if (isascii) PetscCall(PetscLimiterView_VanLeer_Ascii(lim, viewer));
    642:   PetscFunctionReturn(PETSC_SUCCESS);
    643: }

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

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

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

    664:   Level: intermediate

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

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

    673:   PetscFunctionBegin;
    675:   PetscCall(PetscNew(&l));
    676:   lim->data = l;

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

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

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

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

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

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

    705:   PetscFunctionBegin;
    708:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    709:   if (isascii) PetscCall(PetscLimiterView_VanAlbada_Ascii(lim, viewer));
    710:   PetscFunctionReturn(PETSC_SUCCESS);
    711: }

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

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

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

    732:   Level: intermediate

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

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

    741:   PetscFunctionBegin;
    743:   PetscCall(PetscNew(&l));
    744:   lim->data = l;

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

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

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

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

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

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

    773:   PetscFunctionBegin;
    776:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    777:   if (isascii) PetscCall(PetscLimiterView_Superbee_Ascii(lim, viewer));
    778:   PetscFunctionReturn(PETSC_SUCCESS);
    779: }

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

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

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

    800:   Level: intermediate

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

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

    809:   PetscFunctionBegin;
    811:   PetscCall(PetscNew(&l));
    812:   lim->data = l;

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

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

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

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

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

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

    841:   PetscFunctionBegin;
    844:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    845:   if (isascii) PetscCall(PetscLimiterView_MC_Ascii(lim, viewer));
    846:   PetscFunctionReturn(PETSC_SUCCESS);
    847: }

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

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

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

    869:   Level: intermediate

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

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

    878:   PetscFunctionBegin;
    880:   PetscCall(PetscNew(&l));
    881:   lim->data = l;

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

    887: PetscClassId PETSCFV_CLASSID = 0;

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

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

    895:   Not Collective, No Fortran Support

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

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

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

    916:   Level: advanced

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

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

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

    933:   Collective

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

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

    942:   Level: intermediate

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

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

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

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

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

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

    971:   Not Collective

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

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

    979:   Level: intermediate

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

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

    996:   Collective

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

    1003:   Level: intermediate

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

    1015: /*@
    1016:   PetscFVView - Views a `PetscFV`

    1018:   Collective

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

    1024:   Level: beginner

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

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

    1040:   Collective

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

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

    1048:   Level: intermediate

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

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

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

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

    1084:   Collective

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

    1089:   Level: intermediate

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

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

    1105:   Collective

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

    1110:   Level: beginner

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

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

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

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

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

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

    1144:   Collective

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

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

    1152:   Level: beginner

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

    1160:   PetscFunctionBegin;
    1161:   PetscAssertPointer(fvm, 2);
    1162:   PetscCall(PetscFVInitializePackage());

    1164:   PetscCall(PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView));
    1165:   PetscCall(PetscMemzero(f->ops, sizeof(struct _PetscFVOps)));
    1166:   PetscCall(PetscLimiterCreate(comm, &f->limiter));
    1167:   f->numComponents    = 1;
    1168:   f->dim              = 0;
    1169:   f->computeGradients = PETSC_FALSE;
    1170:   f->fluxWork         = NULL;
    1171:   PetscCall(PetscCalloc1(f->numComponents, &f->componentNames));

    1173:   *fvm = f;
    1174:   PetscFunctionReturn(PETSC_SUCCESS);
    1175: }

    1177: /*@
    1178:   PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`

    1180:   Logically Collective

    1182:   Input Parameters:
    1183: + fvm - the `PetscFV` object
    1184: - lim - The `PetscLimiter`

    1186:   Level: intermediate

    1188: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
    1189: @*/
    1190: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
    1191: {
    1192:   PetscFunctionBegin;
    1195:   PetscCall(PetscLimiterDestroy(&fvm->limiter));
    1196:   PetscCall(PetscObjectReference((PetscObject)lim));
    1197:   fvm->limiter = lim;
    1198:   PetscFunctionReturn(PETSC_SUCCESS);
    1199: }

    1201: /*@
    1202:   PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`

    1204:   Not Collective

    1206:   Input Parameter:
    1207: . fvm - the `PetscFV` object

    1209:   Output Parameter:
    1210: . lim - The `PetscLimiter`

    1212:   Level: intermediate

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

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

    1228:   Logically Collective

    1230:   Input Parameters:
    1231: + fvm  - the `PetscFV` object
    1232: - comp - The number of components

    1234:   Level: intermediate

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

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

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

    1258:   Not Collective

    1260:   Input Parameter:
    1261: . fvm - the `PetscFV` object

    1263:   Output Parameter:
    1264: . comp - The number of components

    1266:   Level: intermediate

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

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

    1282:   Logically Collective

    1284:   Input Parameters:
    1285: + fvm  - the `PetscFV` object
    1286: . comp - the component number
    1287: - name - the component name

    1289:   Level: intermediate

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

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

    1304:   Logically Collective

    1306:   Input Parameters:
    1307: + fvm  - the `PetscFV` object
    1308: - comp - the component number

    1310:   Output Parameter:
    1311: . name - the component name

    1313:   Level: intermediate

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

    1324: /*@
    1325:   PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`

    1327:   Logically Collective

    1329:   Input Parameters:
    1330: + fvm - the `PetscFV` object
    1331: - dim - The spatial dimension

    1333:   Level: intermediate

    1335: .seealso: `PetscFV`, `PetscFVGetSpatialDimension()`
    1336: @*/
    1337: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
    1338: {
    1339:   PetscFunctionBegin;
    1341:   fvm->dim = dim;
    1342:   PetscFunctionReturn(PETSC_SUCCESS);
    1343: }

    1345: /*@
    1346:   PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`

    1348:   Not Collective

    1350:   Input Parameter:
    1351: . fvm - the `PetscFV` object

    1353:   Output Parameter:
    1354: . dim - The spatial dimension

    1356:   Level: intermediate

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

    1369: /*@
    1370:   PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`

    1372:   Logically Collective

    1374:   Input Parameters:
    1375: + fvm              - the `PetscFV` object
    1376: - computeGradients - Flag to compute cell gradients

    1378:   Level: intermediate

    1380: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
    1381: @*/
    1382: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
    1383: {
    1384:   PetscFunctionBegin;
    1386:   fvm->computeGradients = computeGradients;
    1387:   PetscFunctionReturn(PETSC_SUCCESS);
    1388: }

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

    1393:   Not Collective

    1395:   Input Parameter:
    1396: . fvm - the `PetscFV` object

    1398:   Output Parameter:
    1399: . computeGradients - Flag to compute cell gradients

    1401:   Level: intermediate

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

    1414: /*@
    1415:   PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`

    1417:   Logically Collective

    1419:   Input Parameters:
    1420: + fvm - the `PetscFV` object
    1421: - q   - The `PetscQuadrature`

    1423:   Level: intermediate

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

    1437: /*@
    1438:   PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`

    1440:   Not Collective

    1442:   Input Parameter:
    1443: . fvm - the `PetscFV` object

    1445:   Output Parameter:
    1446: . q - The `PetscQuadrature`

    1448:   Level: intermediate

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

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

    1471: /*@
    1472:   PetscFVCreateDualSpace - Creates a `PetscDualSpace` appropriate for the `PetscFV`

    1474:   Not Collective

    1476:   Input Parameters:
    1477: + fvm - The `PetscFV` object
    1478: - ct  - The `DMPolytopeType` for the cell

    1480:   Level: intermediate

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

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

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

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

    1519:   Not Collective

    1521:   Input Parameter:
    1522: . fvm - The `PetscFV` object

    1524:   Output Parameter:
    1525: . sp - The `PetscDualSpace` object

    1527:   Level: intermediate

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

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

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

    1549: /*@
    1550:   PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product

    1552:   Not Collective

    1554:   Input Parameters:
    1555: + fvm - The `PetscFV` object
    1556: - sp  - The `PetscDualSpace` object

    1558:   Level: intermediate

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

    1563: .seealso: `PetscFVGetDualSpace()`, `PetscFVCreateDualSpace()`, `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
    1564: @*/
    1565: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
    1566: {
    1567:   PetscFunctionBegin;
    1570:   PetscCall(PetscDualSpaceDestroy(&fvm->dualSpace));
    1571:   fvm->dualSpace = sp;
    1572:   PetscCall(PetscObjectReference((PetscObject)fvm->dualSpace));
    1573:   PetscFunctionReturn(PETSC_SUCCESS);
    1574: }

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

    1579:   Not Collective

    1581:   Input Parameter:
    1582: . fvm - The `PetscFV` object

    1584:   Output Parameter:
    1585: . T - The basis function values and derivatives at quadrature points

    1587:   Level: intermediate

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

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

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

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

    1615:   Not Collective

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

    1624:   Output Parameter:
    1625: . T - The basis function values and derivative at tabulation points

    1627:   Level: intermediate

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

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

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

    1685: /*@
    1686:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

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

    1693:   Output Parameter:
    1694: . grad - the gradient

    1696:   Level: advanced

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

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

    1711:   Not Collective

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

    1723:   Output Parameters:
    1724: + fluxL - the left fluxes for each face
    1725: - fluxR - the right fluxes for each face

    1727:   Level: developer

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

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

    1742:   Input Parameter:
    1743: . fv - The initial `PetscFV`

    1745:   Output Parameter:
    1746: . fvNew - A clone of the `PetscFV`

    1748:   Level: advanced

    1750:   Notes:
    1751:   This is typically used to change the number of components.

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

    1762:   PetscFunctionBegin;
    1763:   PetscCall(PetscFVGetDualSpace(fv, &Q));
    1764:   PetscCall(PetscFVGetQuadrature(fv, &q));
    1765:   PetscCall(PetscDualSpaceGetDM(Q, &K));

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

    1777: /*@
    1778:   PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into
    1779:   smaller copies.

    1781:   Input Parameter:
    1782: . fv - The initial `PetscFV`

    1784:   Output Parameter:
    1785: . fvRef - The refined `PetscFV`

    1787:   Level: advanced

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

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

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

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

    1854: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
    1855: {
    1856:   PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;

    1858:   PetscFunctionBegin;
    1859:   PetscCall(PetscFree(b));
    1860:   PetscFunctionReturn(PETSC_SUCCESS);
    1861: }

    1863: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
    1864: {
    1865:   PetscInt          Nc, c;
    1866:   PetscViewerFormat format;

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

    1879: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
    1880: {
    1881:   PetscBool isascii;

    1883:   PetscFunctionBegin;
    1886:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1887:   if (isascii) PetscCall(PetscFVView_Upwind_Ascii(fv, viewer));
    1888:   PetscFunctionReturn(PETSC_SUCCESS);
    1889: }

    1891: static PetscErrorCode PetscFVComputeGradient_Upwind(PetscFV fv, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
    1892: {
    1893:   PetscInt dim;

    1895:   PetscFunctionBegin;
    1896:   PetscCall(PetscFVGetSpatialDimension(fv, &dim));
    1897:   for (PetscInt f = 0; f < numFaces; ++f) {
    1898:     for (PetscInt d = 0; d < dim; ++d) grad[f * dim + d] = 0.;
    1899:   }
    1900:   PetscFunctionReturn(PETSC_SUCCESS);
    1901: }

    1903: /*
    1904:   neighborVol[f*2+0] contains the left  geom
    1905:   neighborVol[f*2+1] contains the right geom
    1906: */
    1907: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
    1908: {
    1909:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
    1910:   void              *rctx;
    1911:   PetscScalar       *flux = fvm->fluxWork;
    1912:   const PetscScalar *constants;
    1913:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;

    1915:   PetscFunctionBegin;
    1916:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
    1917:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
    1918:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
    1919:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
    1920:   PetscCall(PetscDSGetContext(prob, field, &rctx));
    1921:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
    1922:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
    1923:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
    1924:   for (f = 0; f < Nf; ++f) {
    1925:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
    1926:     for (d = 0; d < pdim; ++d) {
    1927:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
    1928:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
    1929:     }
    1930:   }
    1931:   PetscFunctionReturn(PETSC_SUCCESS);
    1932: }

    1934: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
    1935: {
    1936:   PetscFunctionBegin;
    1937:   fvm->ops->setfromoptions       = NULL;
    1938:   fvm->ops->view                 = PetscFVView_Upwind;
    1939:   fvm->ops->destroy              = PetscFVDestroy_Upwind;
    1940:   fvm->ops->computegradient      = PetscFVComputeGradient_Upwind;
    1941:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
    1942:   PetscFunctionReturn(PETSC_SUCCESS);
    1943: }

    1945: /*MC
    1946:   PETSCFVUPWIND = "upwind" - A `PetscFV` implementation

    1948:   Level: intermediate

    1950: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
    1951: M*/

    1953: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
    1954: {
    1955:   PetscFV_Upwind *b;

    1957:   PetscFunctionBegin;
    1959:   PetscCall(PetscNew(&b));
    1960:   fvm->data = b;

    1962:   PetscCall(PetscFVInitialize_Upwind(fvm));
    1963:   PetscFunctionReturn(PETSC_SUCCESS);
    1964: }

    1966: #include <petscblaslapack.h>

    1968: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
    1969: {
    1970:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;

    1972:   PetscFunctionBegin;
    1973:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL));
    1974:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
    1975:   PetscCall(PetscFree(ls));
    1976:   PetscFunctionReturn(PETSC_SUCCESS);
    1977: }

    1979: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
    1980: {
    1981:   PetscInt          Nc, c;
    1982:   PetscViewerFormat format;

    1984:   PetscFunctionBegin;
    1985:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
    1986:   PetscCall(PetscViewerGetFormat(viewer, &format));
    1987:   PetscCall(PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n"));
    1988:   PetscCall(PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc));
    1989:   for (c = 0; c < Nc; c++) {
    1990:     if (fv->componentNames[c]) PetscCall(PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]));
    1991:   }
    1992:   PetscFunctionReturn(PETSC_SUCCESS);
    1993: }

    1995: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
    1996: {
    1997:   PetscBool isascii;

    1999:   PetscFunctionBegin;
    2002:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    2003:   if (isascii) PetscCall(PetscFVView_LeastSquares_Ascii(fv, viewer));
    2004:   PetscFunctionReturn(PETSC_SUCCESS);
    2005: }

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

    2014:   PetscFunctionBegin;
    2015:   if (debug) {
    2016:     PetscCall(PetscMalloc1(m * n, &Aback));
    2017:     PetscCall(PetscArraycpy(Aback, A, m * n));
    2018:   }

    2020:   PetscCall(PetscBLASIntCast(m, &M));
    2021:   PetscCall(PetscBLASIntCast(n, &N));
    2022:   PetscCall(PetscBLASIntCast(mstride, &lda));
    2023:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
    2024:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    2025:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
    2026:   PetscCall(PetscFPTrapPop());
    2027:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGEQRF error");
    2028:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

    2030:   /* Extract an explicit representation of Q */
    2031:   Q = Ainv;
    2032:   PetscCall(PetscArraycpy(Q, A, mstride * n));
    2033:   K = N; /* full rank */
    2034:   PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));
    2035:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xORGQR/xUNGQR error");

    2037:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
    2038:   Alpha = 1.0;
    2039:   ldb   = lda;
    2040:   BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
    2041:   /* Ainv is Q, overwritten with inverse */

    2043:   if (debug) { /* Check that pseudo-inverse worked */
    2044:     PetscScalar  Beta = 0.0;
    2045:     PetscBLASInt ldc;
    2046:     K   = N;
    2047:     ldc = N;
    2048:     BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
    2049:     PetscCall(PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF));
    2050:     PetscCall(PetscFree(Aback));
    2051:   }
    2052:   PetscFunctionReturn(PETSC_SUCCESS);
    2053: }

    2055: /* Overwrites A. Can handle degenerate problems and m<n. */
    2056: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
    2057: {
    2058:   PetscScalar *Brhs;
    2059:   PetscScalar *tmpwork;
    2060:   PetscReal    rcond;
    2061: #if defined(PETSC_USE_COMPLEX)
    2062:   PetscInt   rworkSize;
    2063:   PetscReal *rwork, *rtau;
    2064: #endif
    2065:   PetscInt     i, j, maxmn;
    2066:   PetscBLASInt M, N, lda, ldb, ldwork;
    2067:   PetscBLASInt nrhs, irank, info;

    2069:   PetscFunctionBegin;
    2070:   /* initialize to identity */
    2071:   tmpwork = work;
    2072:   Brhs    = Ainv;
    2073:   maxmn   = PetscMax(m, n);
    2074:   for (j = 0; j < maxmn; j++) {
    2075:     for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
    2076:   }

    2078:   PetscCall(PetscBLASIntCast(m, &M));
    2079:   PetscCall(PetscBLASIntCast(n, &N));
    2080:   PetscCall(PetscBLASIntCast(mstride, &lda));
    2081:   PetscCall(PetscBLASIntCast(maxmn, &ldb));
    2082:   PetscCall(PetscBLASIntCast(worksize, &ldwork));
    2083:   rcond = -1;
    2084:   nrhs  = M;
    2085: #if defined(PETSC_USE_COMPLEX)
    2086:   rworkSize = 5 * PetscMin(M, N);
    2087:   PetscCall(PetscMalloc1(rworkSize, &rwork));
    2088:   PetscCall(PetscMalloc1(PetscMin(M, N), &rtau));
    2089:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    2090:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, rtau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
    2091:   PetscCall(PetscFPTrapPop());
    2092:   PetscCall(PetscFree(rwork));
    2093:   for (i = 0; i < PetscMin(M, N); i++) tau[i] = rtau[i];
    2094:   PetscCall(PetscFree(rtau));
    2095: #else
    2096:   nrhs = M;
    2097:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    2098:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, tau, &rcond, &irank, tmpwork, &ldwork, &info));
    2099:   PetscCall(PetscFPTrapPop());
    2100: #endif
    2101:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "xGELSS error");
    2102:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
    2103:   PetscCheck(irank >= PetscMin(M, N), PETSC_COMM_SELF, PETSC_ERR_USER, "Rank deficient least squares fit, indicates an isolated cell with two collinear points");
    2104:   PetscFunctionReturn(PETSC_SUCCESS);
    2105: }

    2107: #if 0
    2108: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
    2109: {
    2110:   PetscReal       grad[2] = {0, 0};
    2111:   const PetscInt *faces;
    2112:   PetscInt        numFaces, f;

    2114:   PetscFunctionBegin;
    2115:   PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
    2116:   PetscCall(DMPlexGetCone(dm, cell, &faces));
    2117:   for (f = 0; f < numFaces; ++f) {
    2118:     const PetscInt *fcells;
    2119:     const CellGeom *cg1;
    2120:     const FaceGeom *fg;

    2122:     PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
    2123:     PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg));
    2124:     for (i = 0; i < 2; ++i) {
    2125:       PetscScalar du;

    2127:       if (fcells[i] == c) continue;
    2128:       PetscCall(DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1));
    2129:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
    2130:       grad[0] += fg->grad[!i][0] * du;
    2131:       grad[1] += fg->grad[!i][1] * du;
    2132:     }
    2133:   }
    2134:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]));
    2135:   PetscFunctionReturn(PETSC_SUCCESS);
    2136: }
    2137: #endif

    2139: /*
    2140:   PetscFVComputeGradient_LeastSquares - Compute the gradient reconstruction matrix for a given cell

    2142:   Input Parameters:
    2143: + fvm      - The `PetscFV` object
    2144: . numFaces - The number of cell faces which are not constrained
    2145: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

    2147:   Level: developer

    2149: .seealso: `PetscFV`, `PetscFVCreate()`
    2150: */
    2151: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
    2152: {
    2153:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *)fvm->data;
    2154:   const PetscBool       useSVD   = PETSC_TRUE;
    2155:   const PetscInt        maxFaces = ls->maxFaces;
    2156:   PetscInt              dim, f, d;

    2158:   PetscFunctionBegin;
    2159:   if (numFaces > maxFaces) {
    2160:     PetscCheck(maxFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
    2161:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
    2162:   }
    2163:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
    2164:   for (f = 0; f < numFaces; ++f) {
    2165:     for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
    2166:   }
    2167:   /* Overwrites B with garbage, returns Binv in row-major format */
    2168:   if (useSVD) {
    2169:     PetscInt maxmn = PetscMax(numFaces, dim);
    2170:     PetscCall(PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
    2171:     /* Binv shaped in column-major, coldim=maxmn.*/
    2172:     for (f = 0; f < numFaces; ++f) {
    2173:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
    2174:     }
    2175:   } else {
    2176:     PetscCall(PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work));
    2177:     /* Binv shaped in row-major, rowdim=maxFaces.*/
    2178:     for (f = 0; f < numFaces; ++f) {
    2179:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
    2180:     }
    2181:   }
    2182:   PetscFunctionReturn(PETSC_SUCCESS);
    2183: }

    2185: /*
    2186:   neighborVol[f*2+0] contains the left  geom
    2187:   neighborVol[f*2+1] contains the right geom
    2188: */
    2189: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
    2190: {
    2191:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
    2192:   void              *rctx;
    2193:   PetscScalar       *flux = fvm->fluxWork;
    2194:   const PetscScalar *constants;
    2195:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;

    2197:   PetscFunctionBegin;
    2198:   PetscCall(PetscDSGetTotalComponents(prob, &Nc));
    2199:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
    2200:   PetscCall(PetscDSGetFieldOffset(prob, field, &off));
    2201:   PetscCall(PetscDSGetRiemannSolver(prob, field, &riemann));
    2202:   PetscCall(PetscDSGetContext(prob, field, &rctx));
    2203:   PetscCall(PetscDSGetConstants(prob, &numConstants, &constants));
    2204:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
    2205:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
    2206:   for (f = 0; f < Nf; ++f) {
    2207:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
    2208:     for (d = 0; d < pdim; ++d) {
    2209:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
    2210:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
    2211:     }
    2212:   }
    2213:   PetscFunctionReturn(PETSC_SUCCESS);
    2214: }

    2216: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
    2217: {
    2218:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
    2219:   PetscInt              dim, m, n, nrhs, minmn, maxmn;

    2221:   PetscFunctionBegin;
    2223:   PetscCall(PetscFVGetSpatialDimension(fvm, &dim));
    2224:   PetscCall(PetscFree4(ls->B, ls->Binv, ls->tau, ls->work));
    2225:   ls->maxFaces = maxFaces;
    2226:   m            = ls->maxFaces;
    2227:   n            = dim;
    2228:   nrhs         = ls->maxFaces;
    2229:   minmn        = PetscMin(m, n);
    2230:   maxmn        = PetscMax(m, n);
    2231:   ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
    2232:   PetscCall(PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work));
    2233:   PetscFunctionReturn(PETSC_SUCCESS);
    2234: }

    2236: static PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
    2237: {
    2238:   PetscFunctionBegin;
    2239:   fvm->ops->setfromoptions       = NULL;
    2240:   fvm->ops->view                 = PetscFVView_LeastSquares;
    2241:   fvm->ops->destroy              = PetscFVDestroy_LeastSquares;
    2242:   fvm->ops->computegradient      = PetscFVComputeGradient_LeastSquares;
    2243:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
    2244:   PetscFunctionReturn(PETSC_SUCCESS);
    2245: }

    2247: /*MC
    2248:   PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation

    2250:   Level: intermediate

    2252: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
    2253: M*/

    2255: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
    2256: {
    2257:   PetscFV_LeastSquares *ls;

    2259:   PetscFunctionBegin;
    2261:   PetscCall(PetscNew(&ls));
    2262:   fvm->data = ls;

    2264:   ls->maxFaces = -1;
    2265:   ls->workSize = -1;
    2266:   ls->B        = NULL;
    2267:   ls->Binv     = NULL;
    2268:   ls->tau      = NULL;
    2269:   ls->work     = NULL;

    2271:   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
    2272:   PetscCall(PetscFVInitialize_LeastSquares(fvm));
    2273:   PetscCall(PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS));
    2274:   PetscFunctionReturn(PETSC_SUCCESS);
    2275: }

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

    2280:   Not Collective

    2282:   Input Parameters:
    2283: + fvm      - The `PetscFV` object
    2284: - maxFaces - The maximum number of cell faces

    2286:   Level: intermediate

    2288: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
    2289: @*/
    2290: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
    2291: {
    2292:   PetscFunctionBegin;
    2294:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
    2295:   PetscFunctionReturn(PETSC_SUCCESS);
    2296: }