Actual source code: fv.c

  1: #include <petsc/private/petscfvimpl.h>
  2: #include <petscdmplex.h>
  3: #include <petscdmplextransform.h>
  4: #include <petscds.h>

  6: PetscClassId PETSCLIMITER_CLASSID = 0;

  8: PetscFunctionList PetscLimiterList              = NULL;
  9: PetscBool         PetscLimiterRegisterAllCalled = PETSC_FALSE;

 11: PetscBool  Limitercite       = PETSC_FALSE;
 12: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
 13:                                "  title   = {Analysis of slope limiters on irregular grids},\n"
 14:                                "  journal = {AIAA paper},\n"
 15:                                "  author  = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
 16:                                "  volume  = {490},\n"
 17:                                "  year    = {2005}\n}\n";

 19: /*@C
 20:   PetscLimiterRegister - Adds a new `PetscLimiter` implementation

 22:   Not Collective

 24:   Input Parameters:
 25: + name        - The name of a new user-defined creation routine
 26: - create_func - The creation routine itself

 28:   Sample 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:   PetscFunctionListAdd(&PetscLimiterList, sname, function);
 53:   return 0;
 54: }

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

 59:   Collective on lim

 61:   Input Parameters:
 62: + lim  - The `PetscLimiter` object
 63: - name - The kind of limiter

 65:   Options Database Key:
 66: . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types

 68:   Level: intermediate

 70: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterGetType()`, `PetscLimiterCreate()`
 71: @*/
 72: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
 73: {
 74:   PetscErrorCode (*r)(PetscLimiter);
 75:   PetscBool match;

 78:   PetscObjectTypeCompare((PetscObject)lim, name, &match);
 79:   if (match) return 0;

 81:   PetscLimiterRegisterAll();
 82:   PetscFunctionListFind(PetscLimiterList, name, &r);

 85:   if (lim->ops->destroy) {
 86:     PetscUseTypeMethod(lim, destroy);
 87:     lim->ops->destroy = NULL;
 88:   }
 89:   (*r)(lim);
 90:   PetscObjectChangeTypeName((PetscObject)lim, name);
 91:   return 0;
 92: }

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

 97:   Not Collective

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

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

105:   Level: intermediate

107: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
108: @*/
109: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
110: {
113:   PetscLimiterRegisterAll();
114:   *name = ((PetscObject)lim)->type_name;
115:   return 0;
116: }

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

121:    Collective on A

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

128:    Level: intermediate

130: .seealso: `PetscLimiter`, `PetscLimiterView()`, `PetscObjectViewFromOptions()`, `PetscLimiterCreate()`
131: @*/
132: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A, PetscObject obj, const char name[])
133: {
135:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
136:   return 0;
137: }

139: /*@C
140:   PetscLimiterView - Views a `PetscLimiter`

142:   Collective on lim

144:   Input Parameters:
145: + lim - the `PetscLimiter` object to view
146: - v   - the viewer

148:   Level: beginner

150: .seealso: `PetscLimiter`, `PetscViewer`, `PetscLimiterDestroy()`, `PetscLimiterViewFromOptions()`
151: @*/
152: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
153: {
155:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)lim), &v);
156:   PetscTryTypeMethod(lim, view, v);
157:   return 0;
158: }

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

163:   Collective on lim

165:   Input Parameter:
166: . lim - the `PetscLimiter` object to set options for

168:   Level: intermediate

170: .seealso: `PetscLimiter`, ``PetscLimiterView()`
171: @*/
172: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
173: {
174:   const char *defaultType;
175:   char        name[256];
176:   PetscBool   flg;

179:   if (!((PetscObject)lim)->type_name) defaultType = PETSCLIMITERSIN;
180:   else defaultType = ((PetscObject)lim)->type_name;
181:   PetscLimiterRegisterAll();

183:   PetscObjectOptionsBegin((PetscObject)lim);
184:   PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
185:   if (flg) {
186:     PetscLimiterSetType(lim, name);
187:   } else if (!((PetscObject)lim)->type_name) {
188:     PetscLimiterSetType(lim, defaultType);
189:   }
190:   PetscTryTypeMethod(lim, setfromoptions);
191:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
192:   PetscObjectProcessOptionsHandlers((PetscObject)lim, PetscOptionsObject);
193:   PetscOptionsEnd();
194:   PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
195:   return 0;
196: }

198: /*@C
199:   PetscLimiterSetUp - Construct data structures for the `PetscLimiter`

201:   Collective on lim

203:   Input Parameter:
204: . lim - the `PetscLimiter` object to setup

206:   Level: intermediate

208: .seealso: `PetscLimiter`, ``PetscLimiterView()`, `PetscLimiterDestroy()`
209: @*/
210: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
211: {
213:   PetscTryTypeMethod(lim, setup);
214:   return 0;
215: }

217: /*@
218:   PetscLimiterDestroy - Destroys a `PetscLimiter` object

220:   Collective on lim

222:   Input Parameter:
223: . lim - the `PetscLimiter` object to destroy

225:   Level: beginner

227: .seealso: `PetscLimiter`, `PetscLimiterView()`
228: @*/
229: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
230: {
231:   if (!*lim) return 0;

234:   if (--((PetscObject)(*lim))->refct > 0) {
235:     *lim = NULL;
236:     return 0;
237:   }
238:   ((PetscObject)(*lim))->refct = 0;

240:   PetscTryTypeMethod((*lim), destroy);
241:   PetscHeaderDestroy(lim);
242:   return 0;
243: }

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

248:   Collective

250:   Input Parameter:
251: . comm - The communicator for the `PetscLimiter` object

253:   Output Parameter:
254: . lim - The `PetscLimiter` object

256:   Level: beginner

258: .seealso: `PetscLimiter`, PetscLimiterType`, `PetscLimiterSetType()`, `PETSCLIMITERSIN`
259: @*/
260: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
261: {
262:   PetscLimiter l;

265:   PetscCitationsRegister(LimiterCitation, &Limitercite);
266:   *lim = NULL;
267:   PetscFVInitializePackage();

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

271:   *lim = l;
272:   return 0;
273: }

275: /*@
276:   PetscLimiterLimit - Limit the flux

278:   Input Parameters:
279: + lim  - The `PetscLimiter`
280: - flim - The input field

282:   Output Parameter:
283: . phi  - The limited field

285:   Level: beginner

287:   Note:
288:   Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
289: .vb
290:  The classical flux-limited formulation is psi(r) where

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

294:  The second order TVD region is bounded by

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

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

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

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

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

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

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

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

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

315:  becomes

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

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

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

324: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterSetType()`, `PetscLimiterCreate()`
325: @*/
326: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
327: {
330:   PetscUseTypeMethod(lim, limit, flim, phi);
331:   return 0;
332: }

334: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
335: {
336:   PetscLimiter_Sin *l = (PetscLimiter_Sin *)lim->data;

338:   PetscFree(l);
339:   return 0;
340: }

342: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
343: {
344:   PetscViewerFormat format;

346:   PetscViewerGetFormat(viewer, &format);
347:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
348:   return 0;
349: }

351: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
352: {
353:   PetscBool iascii;

357:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
358:   if (iascii) PetscLimiterView_Sin_Ascii(lim, viewer);
359:   return 0;
360: }

362: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
363: {
364:   *phi = PetscSinReal(PETSC_PI * PetscMax(0, PetscMin(f, 1)));
365:   return 0;
366: }

368: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
369: {
370:   lim->ops->view    = PetscLimiterView_Sin;
371:   lim->ops->destroy = PetscLimiterDestroy_Sin;
372:   lim->ops->limit   = PetscLimiterLimit_Sin;
373:   return 0;
374: }

376: /*MC
377:   PETSCLIMITERSIN = "sin" - A `PetscLimiter` implementation

379:   Level: intermediate

381: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
382: M*/

384: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
385: {
386:   PetscLimiter_Sin *l;

389:   PetscNew(&l);
390:   lim->data = l;

392:   PetscLimiterInitialize_Sin(lim);
393:   return 0;
394: }

396: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
397: {
398:   PetscLimiter_Zero *l = (PetscLimiter_Zero *)lim->data;

400:   PetscFree(l);
401:   return 0;
402: }

404: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
405: {
406:   PetscViewerFormat format;

408:   PetscViewerGetFormat(viewer, &format);
409:   PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
410:   return 0;
411: }

413: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
414: {
415:   PetscBool iascii;

419:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
420:   if (iascii) PetscLimiterView_Zero_Ascii(lim, viewer);
421:   return 0;
422: }

424: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
425: {
426:   *phi = 0.0;
427:   return 0;
428: }

430: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
431: {
432:   lim->ops->view    = PetscLimiterView_Zero;
433:   lim->ops->destroy = PetscLimiterDestroy_Zero;
434:   lim->ops->limit   = PetscLimiterLimit_Zero;
435:   return 0;
436: }

438: /*MC
439:   PETSCLIMITERZERO = "zero" - A simple `PetscLimiter` implementation

441:   Level: intermediate

443: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
444: M*/

446: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
447: {
448:   PetscLimiter_Zero *l;

451:   PetscNew(&l);
452:   lim->data = l;

454:   PetscLimiterInitialize_Zero(lim);
455:   return 0;
456: }

458: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
459: {
460:   PetscLimiter_None *l = (PetscLimiter_None *)lim->data;

462:   PetscFree(l);
463:   return 0;
464: }

466: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
467: {
468:   PetscViewerFormat format;

470:   PetscViewerGetFormat(viewer, &format);
471:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
472:   return 0;
473: }

475: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
476: {
477:   PetscBool iascii;

481:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
482:   if (iascii) PetscLimiterView_None_Ascii(lim, viewer);
483:   return 0;
484: }

486: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
487: {
488:   *phi = 1.0;
489:   return 0;
490: }

492: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
493: {
494:   lim->ops->view    = PetscLimiterView_None;
495:   lim->ops->destroy = PetscLimiterDestroy_None;
496:   lim->ops->limit   = PetscLimiterLimit_None;
497:   return 0;
498: }

500: /*MC
501:   PETSCLIMITERNONE = "none" - A trivial `PetscLimiter` implementation

503:   Level: intermediate

505: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
506: M*/

508: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
509: {
510:   PetscLimiter_None *l;

513:   PetscNew(&l);
514:   lim->data = l;

516:   PetscLimiterInitialize_None(lim);
517:   return 0;
518: }

520: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
521: {
522:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *)lim->data;

524:   PetscFree(l);
525:   return 0;
526: }

528: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
529: {
530:   PetscViewerFormat format;

532:   PetscViewerGetFormat(viewer, &format);
533:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
534:   return 0;
535: }

537: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
538: {
539:   PetscBool iascii;

543:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
544:   if (iascii) PetscLimiterView_Minmod_Ascii(lim, viewer);
545:   return 0;
546: }

548: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
549: {
550:   *phi = 2 * PetscMax(0, PetscMin(f, 1 - f));
551:   return 0;
552: }

554: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
555: {
556:   lim->ops->view    = PetscLimiterView_Minmod;
557:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
558:   lim->ops->limit   = PetscLimiterLimit_Minmod;
559:   return 0;
560: }

562: /*MC
563:   PETSCLIMITERMINMOD = "minmod" - A `PetscLimiter` implementation

565:   Level: intermediate

567: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
568: M*/

570: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
571: {
572:   PetscLimiter_Minmod *l;

575:   PetscNew(&l);
576:   lim->data = l;

578:   PetscLimiterInitialize_Minmod(lim);
579:   return 0;
580: }

582: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
583: {
584:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *)lim->data;

586:   PetscFree(l);
587:   return 0;
588: }

590: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
591: {
592:   PetscViewerFormat format;

594:   PetscViewerGetFormat(viewer, &format);
595:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
596:   return 0;
597: }

599: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
600: {
601:   PetscBool iascii;

605:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
606:   if (iascii) PetscLimiterView_VanLeer_Ascii(lim, viewer);
607:   return 0;
608: }

610: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
611: {
612:   *phi = PetscMax(0, 4 * f * (1 - f));
613:   return 0;
614: }

616: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
617: {
618:   lim->ops->view    = PetscLimiterView_VanLeer;
619:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
620:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
621:   return 0;
622: }

624: /*MC
625:   PETSCLIMITERVANLEER = "vanleer" - A `PetscLimiter` implementation

627:   Level: intermediate

629: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
630: M*/

632: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
633: {
634:   PetscLimiter_VanLeer *l;

637:   PetscNew(&l);
638:   lim->data = l;

640:   PetscLimiterInitialize_VanLeer(lim);
641:   return 0;
642: }

644: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
645: {
646:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *)lim->data;

648:   PetscFree(l);
649:   return 0;
650: }

652: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
653: {
654:   PetscViewerFormat format;

656:   PetscViewerGetFormat(viewer, &format);
657:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
658:   return 0;
659: }

661: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
662: {
663:   PetscBool iascii;

667:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
668:   if (iascii) PetscLimiterView_VanAlbada_Ascii(lim, viewer);
669:   return 0;
670: }

672: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
673: {
674:   *phi = PetscMax(0, 2 * f * (1 - f) / (PetscSqr(f) + PetscSqr(1 - f)));
675:   return 0;
676: }

678: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
679: {
680:   lim->ops->view    = PetscLimiterView_VanAlbada;
681:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
682:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
683:   return 0;
684: }

686: /*MC
687:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter implementation

689:   Level: intermediate

691: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
692: M*/

694: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
695: {
696:   PetscLimiter_VanAlbada *l;

699:   PetscNew(&l);
700:   lim->data = l;

702:   PetscLimiterInitialize_VanAlbada(lim);
703:   return 0;
704: }

706: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
707: {
708:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *)lim->data;

710:   PetscFree(l);
711:   return 0;
712: }

714: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
715: {
716:   PetscViewerFormat format;

718:   PetscViewerGetFormat(viewer, &format);
719:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
720:   return 0;
721: }

723: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
724: {
725:   PetscBool iascii;

729:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
730:   if (iascii) PetscLimiterView_Superbee_Ascii(lim, viewer);
731:   return 0;
732: }

734: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
735: {
736:   *phi = 4 * PetscMax(0, PetscMin(f, 1 - f));
737:   return 0;
738: }

740: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
741: {
742:   lim->ops->view    = PetscLimiterView_Superbee;
743:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
744:   lim->ops->limit   = PetscLimiterLimit_Superbee;
745:   return 0;
746: }

748: /*MC
749:   PETSCLIMITERSUPERBEE = "superbee" - A `PetscLimiter` implementation

751:   Level: intermediate

753: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
754: M*/

756: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
757: {
758:   PetscLimiter_Superbee *l;

761:   PetscNew(&l);
762:   lim->data = l;

764:   PetscLimiterInitialize_Superbee(lim);
765:   return 0;
766: }

768: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
769: {
770:   PetscLimiter_MC *l = (PetscLimiter_MC *)lim->data;

772:   PetscFree(l);
773:   return 0;
774: }

776: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
777: {
778:   PetscViewerFormat format;

780:   PetscViewerGetFormat(viewer, &format);
781:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
782:   return 0;
783: }

785: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
786: {
787:   PetscBool iascii;

791:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
792:   if (iascii) PetscLimiterView_MC_Ascii(lim, viewer);
793:   return 0;
794: }

796: /* aka Barth-Jespersen */
797: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
798: {
799:   *phi = PetscMin(1, 4 * PetscMax(0, PetscMin(f, 1 - f)));
800:   return 0;
801: }

803: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
804: {
805:   lim->ops->view    = PetscLimiterView_MC;
806:   lim->ops->destroy = PetscLimiterDestroy_MC;
807:   lim->ops->limit   = PetscLimiterLimit_MC;
808:   return 0;
809: }

811: /*MC
812:   PETSCLIMITERMC = "mc" - A `PetscLimiter` implementation

814:   Level: intermediate

816: .seealso: `PetscLimiter`, `PetscLimiterType`, `PetscLimiterCreate()`, `PetscLimiterSetType()`
817: M*/

819: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
820: {
821:   PetscLimiter_MC *l;

824:   PetscNew(&l);
825:   lim->data = l;

827:   PetscLimiterInitialize_MC(lim);
828:   return 0;
829: }

831: PetscClassId PETSCFV_CLASSID = 0;

833: PetscFunctionList PetscFVList              = NULL;
834: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

836: /*@C
837:   PetscFVRegister - Adds a new `PetscFV` implementation

839:   Not Collective

841:   Input Parameters:
842: + name        - The name of a new user-defined creation routine
843: - create_func - The creation routine itself

845:   Sample usage:
846: .vb
847:     PetscFVRegister("my_fv", MyPetscFVCreate);
848: .ve

850:   Then, your PetscFV type can be chosen with the procedural interface via
851: .vb
852:     PetscFVCreate(MPI_Comm, PetscFV *);
853:     PetscFVSetType(PetscFV, "my_fv");
854: .ve
855:    or at runtime via the option
856: .vb
857:     -petscfv_type my_fv
858: .ve

860:   Level: advanced

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

865: .seealso: `PetscFV`, `PetscFVType`, `PetscFVRegisterAll()`, `PetscFVRegisterDestroy()`
866: @*/
867: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
868: {
869:   PetscFunctionListAdd(&PetscFVList, sname, function);
870:   return 0;
871: }

873: /*@C
874:   PetscFVSetType - Builds a particular `PetscFV`

876:   Collective on fvm

878:   Input Parameters:
879: + fvm  - The `PetscFV` object
880: - name - The type of FVM space

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

885:   Level: intermediate

887: .seealso: `PetscFV`, `PetscFVType`, `PetscFVGetType()`, `PetscFVCreate()`
888: @*/
889: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
890: {
891:   PetscErrorCode (*r)(PetscFV);
892:   PetscBool match;

895:   PetscObjectTypeCompare((PetscObject)fvm, name, &match);
896:   if (match) return 0;

898:   PetscFVRegisterAll();
899:   PetscFunctionListFind(PetscFVList, name, &r);

902:   PetscTryTypeMethod(fvm, destroy);
903:   fvm->ops->destroy = NULL;

905:   (*r)(fvm);
906:   PetscObjectChangeTypeName((PetscObject)fvm, name);
907:   return 0;
908: }

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

913:   Not Collective

915:   Input Parameter:
916: . fvm  - The `PetscFV`

918:   Output Parameter:
919: . name - The `PetscFVType` name

921:   Level: intermediate

923: .seealso: `PetscFV`, `PetscFVType`, `PetscFVSetType()`, `PetscFVCreate()`
924: @*/
925: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
926: {
929:   PetscFVRegisterAll();
930:   *name = ((PetscObject)fvm)->type_name;
931:   return 0;
932: }

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

937:    Collective on A

939:    Input Parameters:
940: +  A - the `PetscFV` object
941: .  obj - Optional object that provides the options prefix
942: -  name - command line option name

944:    Level: intermediate

946: .seealso: `PetscFV`, `PetscFVView()`, `PetscObjectViewFromOptions()`, `PetscFVCreate()`
947: @*/
948: PetscErrorCode PetscFVViewFromOptions(PetscFV A, PetscObject obj, const char name[])
949: {
951:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
952:   return 0;
953: }

955: /*@C
956:   PetscFVView - Views a `PetscFV`

958:   Collective on fvm

960:   Input Parameters:
961: + fvm - the `PetscFV` object to view
962: - v   - the viewer

964:   Level: beginner

966: .seealso: `PetscFV`, `PetscViewer`, `PetscFVDestroy()`
967: @*/
968: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
969: {
971:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)fvm), &v);
972:   PetscTryTypeMethod(fvm, view, v);
973:   return 0;
974: }

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

979:   Collective on fvm

981:   Input Parameter:
982: . fvm - the `PetscFV` object to set options for

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

987:   Level: intermediate

989: .seealso: `PetscFV`, `PetscFVView()`
990: @*/
991: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
992: {
993:   const char *defaultType;
994:   char        name[256];
995:   PetscBool   flg;

998:   if (!((PetscObject)fvm)->type_name) defaultType = PETSCFVUPWIND;
999:   else defaultType = ((PetscObject)fvm)->type_name;
1000:   PetscFVRegisterAll();

1002:   PetscObjectOptionsBegin((PetscObject)fvm);
1003:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1004:   if (flg) {
1005:     PetscFVSetType(fvm, name);
1006:   } else if (!((PetscObject)fvm)->type_name) {
1007:     PetscFVSetType(fvm, defaultType);
1008:   }
1009:   PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1010:   PetscTryTypeMethod(fvm, setfromoptions);
1011:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1012:   PetscObjectProcessOptionsHandlers((PetscObject)fvm, PetscOptionsObject);
1013:   PetscLimiterSetFromOptions(fvm->limiter);
1014:   PetscOptionsEnd();
1015:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1016:   return 0;
1017: }

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

1022:   Collective on fvm

1024:   Input Parameter:
1025: . fvm - the `PetscFV` object to setup

1027:   Level: intermediate

1029: .seealso: `PetscFV`, `PetscFVView()`, `PetscFVDestroy()`
1030: @*/
1031: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1032: {
1034:   PetscLimiterSetUp(fvm->limiter);
1035:   PetscTryTypeMethod(fvm, setup);
1036:   return 0;
1037: }

1039: /*@
1040:   PetscFVDestroy - Destroys a `PetscFV` object

1042:   Collective on fvm

1044:   Input Parameter:
1045: . fvm - the `PetscFV` object to destroy

1047:   Level: beginner

1049: .seealso: `PetscFV`, `PetscFVCreate()`, `PetscFVView()`
1050: @*/
1051: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1052: {
1053:   PetscInt i;

1055:   if (!*fvm) return 0;

1058:   if (--((PetscObject)(*fvm))->refct > 0) {
1059:     *fvm = NULL;
1060:     return 0;
1061:   }
1062:   ((PetscObject)(*fvm))->refct = 0;

1064:   for (i = 0; i < (*fvm)->numComponents; i++) PetscFree((*fvm)->componentNames[i]);
1065:   PetscFree((*fvm)->componentNames);
1066:   PetscLimiterDestroy(&(*fvm)->limiter);
1067:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1068:   PetscFree((*fvm)->fluxWork);
1069:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1070:   PetscTabulationDestroy(&(*fvm)->T);

1072:   PetscTryTypeMethod((*fvm), destroy);
1073:   PetscHeaderDestroy(fvm);
1074:   return 0;
1075: }

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

1080:   Collective

1082:   Input Parameter:
1083: . comm - The communicator for the `PetscFV` object

1085:   Output Parameter:
1086: . fvm - The `PetscFV` object

1088:   Level: beginner

1090: .seealso: `PetscFVSet`, `PetscFVSetType()`, `PETSCFVUPWIND`, `PetscFVDestroy()`
1091: @*/
1092: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1093: {
1094:   PetscFV f;

1097:   *fvm = NULL;
1098:   PetscFVInitializePackage();

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

1103:   PetscLimiterCreate(comm, &f->limiter);
1104:   f->numComponents    = 1;
1105:   f->dim              = 0;
1106:   f->computeGradients = PETSC_FALSE;
1107:   f->fluxWork         = NULL;
1108:   PetscCalloc1(f->numComponents, &f->componentNames);

1110:   *fvm = f;
1111:   return 0;
1112: }

1114: /*@
1115:   PetscFVSetLimiter - Set the `PetscLimiter` to the `PetscFV`

1117:   Logically collective on fvm

1119:   Input Parameters:
1120: + fvm - the `PetscFV` object
1121: - lim - The `PetscLimiter`

1123:   Level: intermediate

1125: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVGetLimiter()`
1126: @*/
1127: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1128: {
1131:   PetscLimiterDestroy(&fvm->limiter);
1132:   PetscObjectReference((PetscObject)lim);
1133:   fvm->limiter = lim;
1134:   return 0;
1135: }

1137: /*@
1138:   PetscFVGetLimiter - Get the `PetscLimiter` object from the `PetscFV`

1140:   Not collective

1142:   Input Parameter:
1143: . fvm - the `PetscFV` object

1145:   Output Parameter:
1146: . lim - The `PetscLimiter`

1148:   Level: intermediate

1150: .seealso: `PetscFV`, `PetscLimiter`, `PetscFVSetLimiter()`
1151: @*/
1152: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1153: {
1156:   *lim = fvm->limiter;
1157:   return 0;
1158: }

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

1163:   Logically collective on fvm

1165:   Input Parameters:
1166: + fvm - the `PetscFV` object
1167: - comp - The number of components

1169:   Level: intermediate

1171: .seealso: `PetscFV`, `PetscFVGetNumComponents()`
1172: @*/
1173: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1174: {
1176:   if (fvm->numComponents != comp) {
1177:     PetscInt i;

1179:     for (i = 0; i < fvm->numComponents; i++) PetscFree(fvm->componentNames[i]);
1180:     PetscFree(fvm->componentNames);
1181:     PetscCalloc1(comp, &fvm->componentNames);
1182:   }
1183:   fvm->numComponents = comp;
1184:   PetscFree(fvm->fluxWork);
1185:   PetscMalloc1(comp, &fvm->fluxWork);
1186:   return 0;
1187: }

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

1192:   Not collective

1194:   Input Parameter:
1195: . fvm - the `PetscFV` object

1197:   Output Parameter:
1198: , comp - The number of components

1200:   Level: intermediate

1202: .seealso: `PetscFV`, `PetscFVSetNumComponents()`, `PetscFVSetComponentName()`
1203: @*/
1204: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1205: {
1208:   *comp = fvm->numComponents;
1209:   return 0;
1210: }

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

1215:   Logically collective on fvm

1217:   Input Parameters:
1218: + fvm - the `PetscFV` object
1219: . comp - the component number
1220: - name - the component name

1222:   Level: intermediate

1224: .seealso: `PetscFV`, `PetscFVGetComponentName()`
1225: @*/
1226: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1227: {
1228:   PetscFree(fvm->componentNames[comp]);
1229:   PetscStrallocpy(name, &fvm->componentNames[comp]);
1230:   return 0;
1231: }

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

1236:   Logically collective on fvm
1237:   Input Parameters:
1238: + fvm - the `PetscFV` object
1239: - comp - the component number

1241:   Output Parameter:
1242: . name - the component name

1244:   Level: intermediate

1246: .seealso: `PetscFV`, `PetscFVSetComponentName()`
1247: @*/
1248: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1249: {
1250:   *name = fvm->componentNames[comp];
1251:   return 0;
1252: }

1254: /*@
1255:   PetscFVSetSpatialDimension - Set the spatial dimension of a `PetscFV`

1257:   Logically collective on fvm

1259:   Input Parameters:
1260: + fvm - the `PetscFV` object
1261: - dim - The spatial dimension

1263:   Level: intermediate

1265: .seealso: `PetscFV`, ``PetscFVGetSpatialDimension()`
1266: @*/
1267: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1268: {
1270:   fvm->dim = dim;
1271:   return 0;
1272: }

1274: /*@
1275:   PetscFVGetSpatialDimension - Get the spatial dimension of a `PetscFV`

1277:   Logically collective on fvm

1279:   Input Parameter:
1280: . fvm - the `PetscFV` object

1282:   Output Parameter:
1283: . dim - The spatial dimension

1285:   Level: intermediate

1287: .seealso: `PetscFV`, `PetscFVSetSpatialDimension()`
1288: @*/
1289: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1290: {
1293:   *dim = fvm->dim;
1294:   return 0;
1295: }

1297: /*@
1298:  PetscFVSetComputeGradients - Toggle computation of cell gradients on a `PetscFV`

1300:   Logically collective on fvm

1302:   Input Parameters:
1303: + fvm - the `PetscFV` object
1304: - computeGradients - Flag to compute cell gradients

1306:   Level: intermediate

1308: .seealso: `PetscFV`, `PetscFVGetComputeGradients()`
1309: @*/
1310: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1311: {
1313:   fvm->computeGradients = computeGradients;
1314:   return 0;
1315: }

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

1320:   Not collective

1322:   Input Parameter:
1323: . fvm - the `PetscFV` object

1325:   Output Parameter:
1326: . computeGradients - Flag to compute cell gradients

1328:   Level: intermediate

1330: .seealso: `PetscFV`, `PetscFVSetComputeGradients()`
1331: @*/
1332: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1333: {
1336:   *computeGradients = fvm->computeGradients;
1337:   return 0;
1338: }

1340: /*@
1341:   PetscFVSetQuadrature - Set the `PetscQuadrature` object for a `PetscFV`

1343:   Logically collective on fvm

1345:   Input Parameters:
1346: + fvm - the `PetscFV` object
1347: - q - The `PetscQuadrature`

1349:   Level: intermediate

1351: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVGetQuadrature()`
1352: @*/
1353: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1354: {
1356:   PetscQuadratureDestroy(&fvm->quadrature);
1357:   PetscObjectReference((PetscObject)q);
1358:   fvm->quadrature = q;
1359:   return 0;
1360: }

1362: /*@
1363:   PetscFVGetQuadrature - Get the `PetscQuadrature` from a `PetscFV`

1365:   Not collective

1367:   Input Parameter:
1368: . fvm - the `PetscFV` object

1370:   Output Parameter:
1371: . lim - The `PetscQuadrature`

1373:   Level: intermediate

1375: .seealso: `PetscQuadrature`, `PetscFV`, `PetscFVSetQuadrature()`
1376: @*/
1377: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1378: {
1381:   if (!fvm->quadrature) {
1382:     /* Create default 1-point quadrature */
1383:     PetscReal *points, *weights;

1385:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1386:     PetscCalloc1(fvm->dim, &points);
1387:     PetscMalloc1(1, &weights);
1388:     weights[0] = 1.0;
1389:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1390:   }
1391:   *q = fvm->quadrature;
1392:   return 0;
1393: }

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

1398:   Not collective

1400:   Input Parameter:
1401: . fvm - The `PetscFV` object

1403:   Output Parameter:
1404: . sp - The PetscDualSpace object

1406:   Level: intermediate

1408:   Developer Note:
1409:   There is overlap between the methods of `PetscFE` and `PetscFV`, they should probably share a common parent class

1411: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1412: @*/
1413: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1414: {
1417:   if (!fvm->dualSpace) {
1418:     DM       K;
1419:     PetscInt dim, Nc, c;

1421:     PetscFVGetSpatialDimension(fvm, &dim);
1422:     PetscFVGetNumComponents(fvm, &Nc);
1423:     PetscDualSpaceCreate(PetscObjectComm((PetscObject)fvm), &fvm->dualSpace);
1424:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1425:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K);
1426:     PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1427:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1428:     DMDestroy(&K);
1429:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1430:     /* Should we be using PetscFVGetQuadrature() here? */
1431:     for (c = 0; c < Nc; ++c) {
1432:       PetscQuadrature qc;
1433:       PetscReal      *points, *weights;

1435:       PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1436:       PetscCalloc1(dim, &points);
1437:       PetscCalloc1(Nc, &weights);
1438:       weights[c] = 1.0;
1439:       PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1440:       PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1441:       PetscQuadratureDestroy(&qc);
1442:     }
1443:     PetscDualSpaceSetUp(fvm->dualSpace);
1444:   }
1445:   *sp = fvm->dualSpace;
1446:   return 0;
1447: }

1449: /*@
1450:   PetscFVSetDualSpace - Sets the `PetscDualSpace` used to define the inner product

1452:   Not collective

1454:   Input Parameters:
1455: + fvm - The `PetscFV` object
1456: - sp  - The `PetscDualSpace` object

1458:   Level: intermediate

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

1463: .seealso: `PetscDualSpace`, `PetscFV`, `PetscFVCreate()`
1464: @*/
1465: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1466: {
1469:   PetscDualSpaceDestroy(&fvm->dualSpace);
1470:   fvm->dualSpace = sp;
1471:   PetscObjectReference((PetscObject)fvm->dualSpace);
1472:   return 0;
1473: }

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

1478:   Not collective

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

1483:   Output Parameter:
1484: . T - The basis function values and derivatives at quadrature points

1486:   Level: intermediate

1488:   Note:
1489: .vb
1490:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1491:   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
1492:   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
1493: .ve

1495: .seealso: `PetscFV`, `PetscTabulation`, `PetscFEGetCellTabulation()`, `PetscFVCreateTabulation()`, `PetscFVGetQuadrature()`, `PetscQuadratureGetData()`
1496: @*/
1497: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1498: {
1499:   PetscInt         npoints;
1500:   const PetscReal *points;

1504:   PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1505:   if (!fvm->T) PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);
1506:   *T = fvm->T;
1507:   return 0;
1508: }

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

1513:   Not collective

1515:   Input Parameters:
1516: + fvm     - The `PetscFV` object
1517: . nrepl   - The number of replicas
1518: . npoints - The number of tabulation points in a replica
1519: . points  - The tabulation point coordinates
1520: - K       - The order of derivative to tabulate

1522:   Output Parameter:
1523: . T - The basis function values and derivative at tabulation points

1525:   Level: intermediate

1527:   Note:
1528: .vb
1529:   T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1530:   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
1531:   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
1532: .ve

1534: .seealso: `PetscFV`, `PetscTabulation`, `PetscFECreateTabulation()`, `PetscTabulationDestroy()`, `PetscFEGetCellTabulation()`
1535: @*/
1536: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1537: {
1538:   PetscInt pdim = 1; /* Dimension of approximation space P */
1539:   PetscInt cdim;     /* Spatial dimension */
1540:   PetscInt Nc;       /* Field components */
1541:   PetscInt k, p, d, c, e;

1543:   if (!npoints || K < 0) {
1544:     *T = NULL;
1545:     return 0;
1546:   }
1550:   PetscFVGetSpatialDimension(fvm, &cdim);
1551:   PetscFVGetNumComponents(fvm, &Nc);
1552:   PetscMalloc1(1, T);
1553:   (*T)->K    = !cdim ? 0 : K;
1554:   (*T)->Nr   = nrepl;
1555:   (*T)->Np   = npoints;
1556:   (*T)->Nb   = pdim;
1557:   (*T)->Nc   = Nc;
1558:   (*T)->cdim = cdim;
1559:   PetscMalloc1((*T)->K + 1, &(*T)->T);
1560:   for (k = 0; k <= (*T)->K; ++k) PetscMalloc1(nrepl * npoints * pdim * Nc * PetscPowInt(cdim, k), &(*T)->T[k]);
1561:   if (K >= 0) {
1562:     for (p = 0; p < nrepl * npoints; ++p)
1563:       for (d = 0; d < pdim; ++d)
1564:         for (c = 0; c < Nc; ++c) (*T)->T[0][(p * pdim + d) * Nc + c] = 1.0;
1565:   }
1566:   if (K >= 1) {
1567:     for (p = 0; p < nrepl * npoints; ++p)
1568:       for (d = 0; d < pdim; ++d)
1569:         for (c = 0; c < Nc; ++c)
1570:           for (e = 0; e < cdim; ++e) (*T)->T[1][((p * pdim + d) * Nc + c) * cdim + e] = 0.0;
1571:   }
1572:   if (K >= 2) {
1573:     for (p = 0; p < nrepl * npoints; ++p)
1574:       for (d = 0; d < pdim; ++d)
1575:         for (c = 0; c < Nc; ++c)
1576:           for (e = 0; e < cdim * cdim; ++e) (*T)->T[2][((p * pdim + d) * Nc + c) * cdim * cdim + e] = 0.0;
1577:   }
1578:   return 0;
1579: }

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

1584:   Input Parameters:
1585: + fvm      - The `PetscFV` object
1586: . numFaces - The number of cell faces which are not constrained
1587: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1589:   Level: advanced

1591: .seealso: `PetscFV`, `PetscFVCreate()`
1592: @*/
1593: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1594: {
1596:   PetscTryTypeMethod(fvm, computegradient, numFaces, dx, grad);
1597:   return 0;
1598: }

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

1603:   Not collective

1605:   Input Parameters:
1606: + fvm          - The `PetscFV` object for the field being integrated
1607: . prob         - The `PetscDS` specifing the discretizations and continuum functions
1608: . field        - The field being integrated
1609: . Nf           - The number of faces in the chunk
1610: . fgeom        - The face geometry for each face in the chunk
1611: . neighborVol  - The volume for each pair of cells in the chunk
1612: . uL           - The state from the cell on the left
1613: - uR           - The state from the cell on the right

1615:   Output Parameters:
1616: + fluxL        - the left fluxes for each face
1617: - fluxR        - the right fluxes for each face

1619:   Level: developer

1621: .seealso: `PetscFV`, `PetscDS`, `PetscFVFaceGeom`, `PetscFVCreate()`
1622: @*/
1623: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1624: {
1626:   PetscTryTypeMethod(fvm, integraterhsfunction, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1627:   return 0;
1628: }

1630: /*@
1631:   PetscFVRefine - Create a "refined" `PetscFV` object that refines the reference cell into smaller copies. This is typically used
1632:   to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1633:   sparsity). It is also used to create an interpolation between regularly refined meshes.

1635:   Input Parameter:
1636: . fv - The initial `PetscFV`

1638:   Output Parameter:
1639: . fvRef - The refined `PetscFV`

1641:   Level: advanced

1643: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1644: @*/
1645: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1646: {
1647:   PetscDualSpace  Q, Qref;
1648:   DM              K, Kref;
1649:   PetscQuadrature q, qref;
1650:   DMPolytopeType  ct;
1651:   DMPlexTransform tr;
1652:   PetscReal      *v0;
1653:   PetscReal      *jac, *invjac;
1654:   PetscInt        numComp, numSubelements, s;

1656:   PetscFVGetDualSpace(fv, &Q);
1657:   PetscFVGetQuadrature(fv, &q);
1658:   PetscDualSpaceGetDM(Q, &K);
1659:   /* Create dual space */
1660:   PetscDualSpaceDuplicate(Q, &Qref);
1661:   DMRefine(K, PetscObjectComm((PetscObject)fv), &Kref);
1662:   PetscDualSpaceSetDM(Qref, Kref);
1663:   DMDestroy(&Kref);
1664:   PetscDualSpaceSetUp(Qref);
1665:   /* Create volume */
1666:   PetscFVCreate(PetscObjectComm((PetscObject)fv), fvRef);
1667:   PetscFVSetDualSpace(*fvRef, Qref);
1668:   PetscFVGetNumComponents(fv, &numComp);
1669:   PetscFVSetNumComponents(*fvRef, numComp);
1670:   PetscFVSetUp(*fvRef);
1671:   /* Create quadrature */
1672:   DMPlexGetCellType(K, 0, &ct);
1673:   DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1674:   DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1675:   DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1676:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1677:   PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1678:   for (s = 0; s < numSubelements; ++s) {
1679:     PetscQuadrature  qs;
1680:     const PetscReal *points, *weights;
1681:     PetscReal       *p, *w;
1682:     PetscInt         dim, Nc, npoints, np;

1684:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1685:     PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1686:     np = npoints / numSubelements;
1687:     PetscMalloc1(np * dim, &p);
1688:     PetscMalloc1(np * Nc, &w);
1689:     PetscArraycpy(p, &points[s * np * dim], np * dim);
1690:     PetscArraycpy(w, &weights[s * np * Nc], np * Nc);
1691:     PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1692:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1693:     PetscQuadratureDestroy(&qs);
1694:   }
1695:   PetscFVSetQuadrature(*fvRef, qref);
1696:   DMPlexTransformDestroy(&tr);
1697:   PetscQuadratureDestroy(&qref);
1698:   PetscDualSpaceDestroy(&Qref);
1699:   return 0;
1700: }

1702: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1703: {
1704:   PetscFV_Upwind *b = (PetscFV_Upwind *)fvm->data;

1706:   PetscFree(b);
1707:   return 0;
1708: }

1710: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1711: {
1712:   PetscInt          Nc, c;
1713:   PetscViewerFormat format;

1715:   PetscFVGetNumComponents(fv, &Nc);
1716:   PetscViewerGetFormat(viewer, &format);
1717:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1718:   PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc);
1719:   for (c = 0; c < Nc; c++) {
1720:     if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1721:   }
1722:   return 0;
1723: }

1725: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1726: {
1727:   PetscBool iascii;

1731:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1732:   if (iascii) PetscFVView_Upwind_Ascii(fv, viewer);
1733:   return 0;
1734: }

1736: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1737: {
1738:   return 0;
1739: }

1741: /*
1742:   neighborVol[f*2+0] contains the left  geom
1743:   neighborVol[f*2+1] contains the right geom
1744: */
1745: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1746: {
1747:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1748:   void              *rctx;
1749:   PetscScalar       *flux = fvm->fluxWork;
1750:   const PetscScalar *constants;
1751:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;

1753:   PetscDSGetTotalComponents(prob, &Nc);
1754:   PetscDSGetTotalDimension(prob, &totDim);
1755:   PetscDSGetFieldOffset(prob, field, &off);
1756:   PetscDSGetRiemannSolver(prob, field, &riemann);
1757:   PetscDSGetContext(prob, field, &rctx);
1758:   PetscDSGetConstants(prob, &numConstants, &constants);
1759:   PetscFVGetSpatialDimension(fvm, &dim);
1760:   PetscFVGetNumComponents(fvm, &pdim);
1761:   for (f = 0; f < Nf; ++f) {
1762:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
1763:     for (d = 0; d < pdim; ++d) {
1764:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
1765:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
1766:     }
1767:   }
1768:   return 0;
1769: }

1771: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1772: {
1773:   fvm->ops->setfromoptions       = NULL;
1774:   fvm->ops->setup                = PetscFVSetUp_Upwind;
1775:   fvm->ops->view                 = PetscFVView_Upwind;
1776:   fvm->ops->destroy              = PetscFVDestroy_Upwind;
1777:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1778:   return 0;
1779: }

1781: /*MC
1782:   PETSCFVUPWIND = "upwind" - A `PetscFV` implementation

1784:   Level: intermediate

1786: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
1787: M*/

1789: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1790: {
1791:   PetscFV_Upwind *b;

1794:   PetscNew(&b);
1795:   fvm->data = b;

1797:   PetscFVInitialize_Upwind(fvm);
1798:   return 0;
1799: }

1801: #include <petscblaslapack.h>

1803: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1804: {
1805:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;

1807:   PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1808:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1809:   PetscFree(ls);
1810:   return 0;
1811: }

1813: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1814: {
1815:   PetscInt          Nc, c;
1816:   PetscViewerFormat format;

1818:   PetscFVGetNumComponents(fv, &Nc);
1819:   PetscViewerGetFormat(viewer, &format);
1820:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1821:   PetscViewerASCIIPrintf(viewer, "  num components: %" PetscInt_FMT "\n", Nc);
1822:   for (c = 0; c < Nc; c++) {
1823:     if (fv->componentNames[c]) PetscViewerASCIIPrintf(viewer, "    component %" PetscInt_FMT ": %s\n", c, fv->componentNames[c]);
1824:   }
1825:   return 0;
1826: }

1828: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1829: {
1830:   PetscBool iascii;

1834:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1835:   if (iascii) PetscFVView_LeastSquares_Ascii(fv, viewer);
1836:   return 0;
1837: }

1839: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1840: {
1841:   return 0;
1842: }

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

1851:   if (debug) {
1852:     PetscMalloc1(m * n, &Aback);
1853:     PetscArraycpy(Aback, A, m * n);
1854:   }

1856:   PetscBLASIntCast(m, &M);
1857:   PetscBLASIntCast(n, &N);
1858:   PetscBLASIntCast(mstride, &lda);
1859:   PetscBLASIntCast(worksize, &ldwork);
1860:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1861:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&M, &N, A, &lda, tau, work, &ldwork, &info));
1862:   PetscFPTrapPop();
1864:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1866:   /* Extract an explicit representation of Q */
1867:   Q = Ainv;
1868:   PetscArraycpy(Q, A, mstride * n);
1869:   K = N; /* full rank */
1870:   PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&M, &N, &K, Q, &lda, tau, work, &ldwork, &info));

1873:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1874:   Alpha = 1.0;
1875:   ldb   = lda;
1876:   BLAStrsm_("Right", "Upper", "ConjugateTranspose", "NotUnitTriangular", &M, &N, &Alpha, R, &lda, Q, &ldb);
1877:   /* Ainv is Q, overwritten with inverse */

1879:   if (debug) { /* Check that pseudo-inverse worked */
1880:     PetscScalar  Beta = 0.0;
1881:     PetscBLASInt ldc;
1882:     K   = N;
1883:     ldc = N;
1884:     BLASgemm_("ConjugateTranspose", "Normal", &N, &K, &M, &Alpha, Ainv, &lda, Aback, &ldb, &Beta, work, &ldc);
1885:     PetscScalarView(n * n, work, PETSC_VIEWER_STDOUT_SELF);
1886:     PetscFree(Aback);
1887:   }
1888:   return 0;
1889: }

1891: /* Overwrites A. Can handle degenerate problems and m<n. */
1892: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m, PetscInt mstride, PetscInt n, PetscScalar *A, PetscScalar *Ainv, PetscScalar *tau, PetscInt worksize, PetscScalar *work)
1893: {
1894:   PetscScalar *Brhs;
1895:   PetscScalar *tmpwork;
1896:   PetscReal    rcond;
1897: #if defined(PETSC_USE_COMPLEX)
1898:   PetscInt   rworkSize;
1899:   PetscReal *rwork;
1900: #endif
1901:   PetscInt     i, j, maxmn;
1902:   PetscBLASInt M, N, lda, ldb, ldwork;
1903:   PetscBLASInt nrhs, irank, info;

1905:   /* initialize to identity */
1906:   tmpwork = work;
1907:   Brhs    = Ainv;
1908:   maxmn   = PetscMax(m, n);
1909:   for (j = 0; j < maxmn; j++) {
1910:     for (i = 0; i < maxmn; i++) Brhs[i + j * maxmn] = 1.0 * (i == j);
1911:   }

1913:   PetscBLASIntCast(m, &M);
1914:   PetscBLASIntCast(n, &N);
1915:   PetscBLASIntCast(mstride, &lda);
1916:   PetscBLASIntCast(maxmn, &ldb);
1917:   PetscBLASIntCast(worksize, &ldwork);
1918:   rcond = -1;
1919:   nrhs  = M;
1920: #if defined(PETSC_USE_COMPLEX)
1921:   rworkSize = 5 * PetscMin(M, N);
1922:   PetscMalloc1(rworkSize, &rwork);
1923:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1924:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, rwork, &info));
1925:   PetscFPTrapPop();
1926:   PetscFree(rwork);
1927: #else
1928:   nrhs = M;
1929:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1930:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &nrhs, A, &lda, Brhs, &ldb, (PetscReal *)tau, &rcond, &irank, tmpwork, &ldwork, &info));
1931:   PetscFPTrapPop();
1932: #endif
1934:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1936:   return 0;
1937: }

1939: #if 0
1940: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1941: {
1942:   PetscReal       grad[2] = {0, 0};
1943:   const PetscInt *faces;
1944:   PetscInt        numFaces, f;

1946:   DMPlexGetConeSize(dm, cell, &numFaces);
1947:   DMPlexGetCone(dm, cell, &faces);
1948:   for (f = 0; f < numFaces; ++f) {
1949:     const PetscInt *fcells;
1950:     const CellGeom *cg1;
1951:     const FaceGeom *fg;

1953:     DMPlexGetSupport(dm, faces[f], &fcells);
1954:     DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1955:     for (i = 0; i < 2; ++i) {
1956:       PetscScalar du;

1958:       if (fcells[i] == c) continue;
1959:       DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1960:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1961:       grad[0] += fg->grad[!i][0] * du;
1962:       grad[1] += fg->grad[!i][1] * du;
1963:     }
1964:   }
1965:   PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1966:   return 0;
1967: }
1968: #endif

1970: /*
1971:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

1973:   Input Parameters:
1974: + fvm      - The `PetscFV` object
1975: . numFaces - The number of cell faces which are not constrained
1976: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1978:   Level: developer

1980: .seealso: `PetscFV`, `PetscFVCreate()`
1981: */
1982: static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1983: {
1984:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *)fvm->data;
1985:   const PetscBool       useSVD   = PETSC_TRUE;
1986:   const PetscInt        maxFaces = ls->maxFaces;
1987:   PetscInt              dim, f, d;

1989:   if (numFaces > maxFaces) {
1991:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %" PetscInt_FMT " > %" PetscInt_FMT " maxfaces", numFaces, maxFaces);
1992:   }
1993:   PetscFVGetSpatialDimension(fvm, &dim);
1994:   for (f = 0; f < numFaces; ++f) {
1995:     for (d = 0; d < dim; ++d) ls->B[d * maxFaces + f] = dx[f * dim + d];
1996:   }
1997:   /* Overwrites B with garbage, returns Binv in row-major format */
1998:   if (useSVD) {
1999:     PetscInt maxmn = PetscMax(numFaces, dim);
2000:     PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2001:     /* Binv shaped in column-major, coldim=maxmn.*/
2002:     for (f = 0; f < numFaces; ++f) {
2003:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d + maxmn * f];
2004:     }
2005:   } else {
2006:     PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
2007:     /* Binv shaped in row-major, rowdim=maxFaces.*/
2008:     for (f = 0; f < numFaces; ++f) {
2009:       for (d = 0; d < dim; ++d) grad[f * dim + d] = ls->Binv[d * maxFaces + f];
2010:     }
2011:   }
2012:   return 0;
2013: }

2015: /*
2016:   neighborVol[f*2+0] contains the left  geom
2017:   neighborVol[f*2+1] contains the right geom
2018: */
2019: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol, PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2020: {
2021:   void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2022:   void              *rctx;
2023:   PetscScalar       *flux = fvm->fluxWork;
2024:   const PetscScalar *constants;
2025:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;

2027:   PetscDSGetTotalComponents(prob, &Nc);
2028:   PetscDSGetTotalDimension(prob, &totDim);
2029:   PetscDSGetFieldOffset(prob, field, &off);
2030:   PetscDSGetRiemannSolver(prob, field, &riemann);
2031:   PetscDSGetContext(prob, field, &rctx);
2032:   PetscDSGetConstants(prob, &numConstants, &constants);
2033:   PetscFVGetSpatialDimension(fvm, &dim);
2034:   PetscFVGetNumComponents(fvm, &pdim);
2035:   for (f = 0; f < Nf; ++f) {
2036:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f * Nc], &uR[f * Nc], numConstants, constants, flux, rctx);
2037:     for (d = 0; d < pdim; ++d) {
2038:       fluxL[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 0];
2039:       fluxR[f * totDim + off + d] = flux[d] / neighborVol[f * 2 + 1];
2040:     }
2041:   }
2042:   return 0;
2043: }

2045: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2046: {
2047:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *)fvm->data;
2048:   PetscInt              dim, m, n, nrhs, minmn, maxmn;

2051:   PetscFVGetSpatialDimension(fvm, &dim);
2052:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2053:   ls->maxFaces = maxFaces;
2054:   m            = ls->maxFaces;
2055:   n            = dim;
2056:   nrhs         = ls->maxFaces;
2057:   minmn        = PetscMin(m, n);
2058:   maxmn        = PetscMax(m, n);
2059:   ls->workSize = 3 * minmn + PetscMax(2 * minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2060:   PetscMalloc4(m * n, &ls->B, maxmn * maxmn, &ls->Binv, minmn, &ls->tau, ls->workSize, &ls->work);
2061:   return 0;
2062: }

2064: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2065: {
2066:   fvm->ops->setfromoptions       = NULL;
2067:   fvm->ops->setup                = PetscFVSetUp_LeastSquares;
2068:   fvm->ops->view                 = PetscFVView_LeastSquares;
2069:   fvm->ops->destroy              = PetscFVDestroy_LeastSquares;
2070:   fvm->ops->computegradient      = PetscFVComputeGradient_LeastSquares;
2071:   fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2072:   return 0;
2073: }

2075: /*MC
2076:   PETSCFVLEASTSQUARES = "leastsquares" - A `PetscFV` implementation

2078:   Level: intermediate

2080: .seealso: `PetscFV`, `PetscFVType`, `PetscFVCreate()`, `PetscFVSetType()`
2081: M*/

2083: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2084: {
2085:   PetscFV_LeastSquares *ls;

2088:   PetscNew(&ls);
2089:   fvm->data = ls;

2091:   ls->maxFaces = -1;
2092:   ls->workSize = -1;
2093:   ls->B        = NULL;
2094:   ls->Binv     = NULL;
2095:   ls->tau      = NULL;
2096:   ls->work     = NULL;

2098:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2099:   PetscFVInitialize_LeastSquares(fvm);
2100:   PetscObjectComposeFunction((PetscObject)fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2101:   return 0;
2102: }

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

2107:   Not collective

2109:   Input parameters:
2110: + fvm      - The `PetscFV` object
2111: - maxFaces - The maximum number of cell faces

2113:   Level: intermediate

2115: .seealso: `PetscFV`, `PetscFVCreate()`, `PETSCFVLEASTSQUARES`, `PetscFVComputeGradient()`
2116: @*/
2117: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2118: {
2120:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV, PetscInt), (fvm, maxFaces));
2121:   return 0;
2122: }