Actual source code: fv.c

petsc-3.12.4 2020-02-04
Report Typos and Errors
  1:  #include <petsc/private/petscfvimpl.h>
  2:  #include <petsc/private/dmpleximpl.h>
  3:  #include <petscds.h>

  5: PetscClassId PETSCLIMITER_CLASSID = 0;

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

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

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

 21:   Not Collective

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

 27:   Notes:
 28:   PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters

 30:   Sample usage:
 31: .vb
 32:     PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
 33: .ve

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

 45:   Level: advanced

 47: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()

 49: @*/
 50: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
 51: {

 55:   PetscFunctionListAdd(&PetscLimiterList, sname, function);
 56:   return(0);
 57: }

 59: /*@C
 60:   PetscLimiterSetType - Builds a particular PetscLimiter

 62:   Collective on lim

 64:   Input Parameters:
 65: + lim  - The PetscLimiter object
 66: - name - The kind of limiter

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

 71:   Level: intermediate

 73: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
 74: @*/
 75: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
 76: {
 77:   PetscErrorCode (*r)(PetscLimiter);
 78:   PetscBool      match;

 83:   PetscObjectTypeCompare((PetscObject) lim, name, &match);
 84:   if (match) return(0);

 86:   PetscLimiterRegisterAll();
 87:   PetscFunctionListFind(PetscLimiterList, name, &r);
 88:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);

 90:   if (lim->ops->destroy) {
 91:     (*lim->ops->destroy)(lim);
 92:     lim->ops->destroy = NULL;
 93:   }
 94:   (*r)(lim);
 95:   PetscObjectChangeTypeName((PetscObject) lim, name);
 96:   return(0);
 97: }

 99: /*@C
100:   PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.

102:   Not Collective

104:   Input Parameter:
105: . lim  - The PetscLimiter

107:   Output Parameter:
108: . name - The PetscLimiter type name

110:   Level: intermediate

112: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
113: @*/
114: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
115: {

121:   PetscLimiterRegisterAll();
122:   *name = ((PetscObject) lim)->type_name;
123:   return(0);
124: }

126: /*@C
127:   PetscLimiterView - Views a PetscLimiter

129:   Collective on lim

131:   Input Parameter:
132: + lim - the PetscLimiter object to view
133: - v   - the viewer

135:   Level: beginner

137: .seealso: PetscLimiterDestroy()
138: @*/
139: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
140: {

145:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
146:   if (lim->ops->view) {(*lim->ops->view)(lim, v);}
147:   return(0);
148: }

150: /*@
151:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

153:   Collective on lim

155:   Input Parameter:
156: . lim - the PetscLimiter object to set options for

158:   Level: intermediate

160: .seealso: PetscLimiterView()
161: @*/
162: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
163: {
164:   const char    *defaultType;
165:   char           name[256];
166:   PetscBool      flg;

171:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
172:   else                                 defaultType = ((PetscObject) lim)->type_name;
173:   PetscLimiterRegisterAll();

175:   PetscObjectOptionsBegin((PetscObject) lim);
176:   PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
177:   if (flg) {
178:     PetscLimiterSetType(lim, name);
179:   } else if (!((PetscObject) lim)->type_name) {
180:     PetscLimiterSetType(lim, defaultType);
181:   }
182:   if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
183:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
184:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
185:   PetscOptionsEnd();
186:   PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
187:   return(0);
188: }

190: /*@C
191:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

193:   Collective on lim

195:   Input Parameter:
196: . lim - the PetscLimiter object to setup

198:   Level: intermediate

200: .seealso: PetscLimiterView(), PetscLimiterDestroy()
201: @*/
202: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
203: {

208:   if (lim->ops->setup) {(*lim->ops->setup)(lim);}
209:   return(0);
210: }

212: /*@
213:   PetscLimiterDestroy - Destroys a PetscLimiter object

215:   Collective on lim

217:   Input Parameter:
218: . lim - the PetscLimiter object to destroy

220:   Level: beginner

222: .seealso: PetscLimiterView()
223: @*/
224: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
225: {

229:   if (!*lim) return(0);

232:   if (--((PetscObject)(*lim))->refct > 0) {*lim = 0; return(0);}
233:   ((PetscObject) (*lim))->refct = 0;

235:   if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
236:   PetscHeaderDestroy(lim);
237:   return(0);
238: }

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

243:   Collective

245:   Input Parameter:
246: . comm - The communicator for the PetscLimiter object

248:   Output Parameter:
249: . lim - The PetscLimiter object

251:   Level: beginner

253: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
254: @*/
255: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
256: {
257:   PetscLimiter   l;

262:   PetscCitationsRegister(LimiterCitation,&Limitercite);
263:   *lim = NULL;
264:   PetscFVInitializePackage();

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

268:   *lim = l;
269:   return(0);
270: }

272: /*@
273:   PetscLimiterLimit - Limit the flux

275:   Input Parameters:
276: + lim  - The PetscLimiter
277: - flim - The input field

279:   Output Parameter:
280: . phi  - The limited field

282: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
283: $ The classical flux-limited formulation is psi(r) where
284: $
285: $ r = (u[0] - u[-1]) / (u[1] - u[0])
286: $
287: $ The second order TVD region is bounded by
288: $
289: $ psi_minmod(r) = min(r,1)      and        psi_superbee(r) = min(2, 2r, max(1,r))
290: $
291: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
292: $ phi(r)(r+1)/2 in which the reconstructed interface values are
293: $
294: $ u(v) = u[0] + phi(r) (grad u)[0] v
295: $
296: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
297: $
298: $ 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))
299: $
300: $ For a nicer symmetric formulation, rewrite in terms of
301: $
302: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
303: $
304: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
305: $
306: $ phi(r) = phi(1/r)
307: $
308: $ becomes
309: $
310: $ w(f) = w(1-f).
311: $
312: $ The limiters below implement this final form w(f). The reference methods are
313: $
314: $ w_minmod(f) = 2 min(f,(1-f))             w_superbee(r) = 4 min((1-f), f)

316:   Level: beginner

318: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
319: @*/
320: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
321: {

327:   (*lim->ops->limit)(lim, flim, phi);
328:   return(0);
329: }

331: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
332: {
333:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
334:   PetscErrorCode    ierr;

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

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

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

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

360:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
361:   if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
362:   return(0);
363: }

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

372: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
373: {
375:   lim->ops->view    = PetscLimiterView_Sin;
376:   lim->ops->destroy = PetscLimiterDestroy_Sin;
377:   lim->ops->limit   = PetscLimiterLimit_Sin;
378:   return(0);
379: }

381: /*MC
382:   PETSCLIMITERSIN = "sin" - A PetscLimiter object

384:   Level: intermediate

386: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
387: M*/

389: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
390: {
391:   PetscLimiter_Sin *l;
392:   PetscErrorCode    ierr;

396:   PetscNewLog(lim, &l);
397:   lim->data = l;

399:   PetscLimiterInitialize_Sin(lim);
400:   return(0);
401: }

403: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
404: {
405:   PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
406:   PetscErrorCode    ierr;

409:   PetscFree(l);
410:   return(0);
411: }

413: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
414: {
415:   PetscViewerFormat format;
416:   PetscErrorCode    ierr;

419:   PetscViewerGetFormat(viewer, &format);
420:   PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
421:   return(0);
422: }

424: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
425: {
426:   PetscBool      iascii;

432:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
433:   if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
434:   return(0);
435: }

437: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
438: {
440:   *phi = 0.0;
441:   return(0);
442: }

444: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
445: {
447:   lim->ops->view    = PetscLimiterView_Zero;
448:   lim->ops->destroy = PetscLimiterDestroy_Zero;
449:   lim->ops->limit   = PetscLimiterLimit_Zero;
450:   return(0);
451: }

453: /*MC
454:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

456:   Level: intermediate

458: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
459: M*/

461: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
462: {
463:   PetscLimiter_Zero *l;
464:   PetscErrorCode     ierr;

468:   PetscNewLog(lim, &l);
469:   lim->data = l;

471:   PetscLimiterInitialize_Zero(lim);
472:   return(0);
473: }

475: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
476: {
477:   PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
478:   PetscErrorCode    ierr;

481:   PetscFree(l);
482:   return(0);
483: }

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

491:   PetscViewerGetFormat(viewer, &format);
492:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
493:   return(0);
494: }

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

504:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
505:   if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
506:   return(0);
507: }

509: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
510: {
512:   *phi = 1.0;
513:   return(0);
514: }

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

525: /*MC
526:   PETSCLIMITERNONE = "none" - A PetscLimiter object

528:   Level: intermediate

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

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

540:   PetscNewLog(lim, &l);
541:   lim->data = l;

543:   PetscLimiterInitialize_None(lim);
544:   return(0);
545: }

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

553:   PetscFree(l);
554:   return(0);
555: }

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

563:   PetscViewerGetFormat(viewer, &format);
564:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
565:   return(0);
566: }

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

576:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
577:   if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
578:   return(0);
579: }

581: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
582: {
584:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
585:   return(0);
586: }

588: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
589: {
591:   lim->ops->view    = PetscLimiterView_Minmod;
592:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
593:   lim->ops->limit   = PetscLimiterLimit_Minmod;
594:   return(0);
595: }

597: /*MC
598:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

600:   Level: intermediate

602: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
603: M*/

605: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
606: {
607:   PetscLimiter_Minmod *l;
608:   PetscErrorCode    ierr;

612:   PetscNewLog(lim, &l);
613:   lim->data = l;

615:   PetscLimiterInitialize_Minmod(lim);
616:   return(0);
617: }

619: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
620: {
621:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
622:   PetscErrorCode    ierr;

625:   PetscFree(l);
626:   return(0);
627: }

629: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
630: {
631:   PetscViewerFormat format;
632:   PetscErrorCode    ierr;

635:   PetscViewerGetFormat(viewer, &format);
636:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
637:   return(0);
638: }

640: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
641: {
642:   PetscBool      iascii;

648:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
649:   if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
650:   return(0);
651: }

653: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
654: {
656:   *phi = PetscMax(0, 4*f*(1-f));
657:   return(0);
658: }

660: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
661: {
663:   lim->ops->view    = PetscLimiterView_VanLeer;
664:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
665:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
666:   return(0);
667: }

669: /*MC
670:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

672:   Level: intermediate

674: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
675: M*/

677: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
678: {
679:   PetscLimiter_VanLeer *l;
680:   PetscErrorCode    ierr;

684:   PetscNewLog(lim, &l);
685:   lim->data = l;

687:   PetscLimiterInitialize_VanLeer(lim);
688:   return(0);
689: }

691: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
692: {
693:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
694:   PetscErrorCode    ierr;

697:   PetscFree(l);
698:   return(0);
699: }

701: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
702: {
703:   PetscViewerFormat format;
704:   PetscErrorCode    ierr;

707:   PetscViewerGetFormat(viewer, &format);
708:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
709:   return(0);
710: }

712: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
713: {
714:   PetscBool      iascii;

720:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
721:   if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
722:   return(0);
723: }

725: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
726: {
728:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
729:   return(0);
730: }

732: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
733: {
735:   lim->ops->view    = PetscLimiterView_VanAlbada;
736:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
737:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
738:   return(0);
739: }

741: /*MC
742:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

744:   Level: intermediate

746: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
747: M*/

749: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
750: {
751:   PetscLimiter_VanAlbada *l;
752:   PetscErrorCode    ierr;

756:   PetscNewLog(lim, &l);
757:   lim->data = l;

759:   PetscLimiterInitialize_VanAlbada(lim);
760:   return(0);
761: }

763: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
764: {
765:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
766:   PetscErrorCode    ierr;

769:   PetscFree(l);
770:   return(0);
771: }

773: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
774: {
775:   PetscViewerFormat format;
776:   PetscErrorCode    ierr;

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

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

792:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
793:   if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
794:   return(0);
795: }

797: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
798: {
800:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
801:   return(0);
802: }

804: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
805: {
807:   lim->ops->view    = PetscLimiterView_Superbee;
808:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
809:   lim->ops->limit   = PetscLimiterLimit_Superbee;
810:   return(0);
811: }

813: /*MC
814:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

816:   Level: intermediate

818: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
819: M*/

821: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
822: {
823:   PetscLimiter_Superbee *l;
824:   PetscErrorCode    ierr;

828:   PetscNewLog(lim, &l);
829:   lim->data = l;

831:   PetscLimiterInitialize_Superbee(lim);
832:   return(0);
833: }

835: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
836: {
837:   PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
838:   PetscErrorCode    ierr;

841:   PetscFree(l);
842:   return(0);
843: }

845: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
846: {
847:   PetscViewerFormat format;
848:   PetscErrorCode    ierr;

851:   PetscViewerGetFormat(viewer, &format);
852:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
853:   return(0);
854: }

856: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
857: {
858:   PetscBool      iascii;

864:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
865:   if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
866:   return(0);
867: }

869: /* aka Barth-Jespersen */
870: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
871: {
873:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
874:   return(0);
875: }

877: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
878: {
880:   lim->ops->view    = PetscLimiterView_MC;
881:   lim->ops->destroy = PetscLimiterDestroy_MC;
882:   lim->ops->limit   = PetscLimiterLimit_MC;
883:   return(0);
884: }

886: /*MC
887:   PETSCLIMITERMC = "mc" - A PetscLimiter object

889:   Level: intermediate

891: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
892: M*/

894: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
895: {
896:   PetscLimiter_MC *l;
897:   PetscErrorCode    ierr;

901:   PetscNewLog(lim, &l);
902:   lim->data = l;

904:   PetscLimiterInitialize_MC(lim);
905:   return(0);
906: }

908: PetscClassId PETSCFV_CLASSID = 0;

910: PetscFunctionList PetscFVList              = NULL;
911: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

913: /*@C
914:   PetscFVRegister - Adds a new PetscFV implementation

916:   Not Collective

918:   Input Parameters:
919: + name        - The name of a new user-defined creation routine
920: - create_func - The creation routine itself

922:   Notes:
923:   PetscFVRegister() may be called multiple times to add several user-defined PetscFVs

925:   Sample usage:
926: .vb
927:     PetscFVRegister("my_fv", MyPetscFVCreate);
928: .ve

930:   Then, your PetscFV type can be chosen with the procedural interface via
931: .vb
932:     PetscFVCreate(MPI_Comm, PetscFV *);
933:     PetscFVSetType(PetscFV, "my_fv");
934: .ve
935:    or at runtime via the option
936: .vb
937:     -petscfv_type my_fv
938: .ve

940:   Level: advanced

942: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()

944: @*/
945: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
946: {

950:   PetscFunctionListAdd(&PetscFVList, sname, function);
951:   return(0);
952: }

954: /*@C
955:   PetscFVSetType - Builds a particular PetscFV

957:   Collective on fvm

959:   Input Parameters:
960: + fvm  - The PetscFV object
961: - name - The kind of FVM space

963:   Options Database Key:
964: . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types

966:   Level: intermediate

968: .seealso: PetscFVGetType(), PetscFVCreate()
969: @*/
970: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
971: {
972:   PetscErrorCode (*r)(PetscFV);
973:   PetscBool      match;

978:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
979:   if (match) return(0);

981:   PetscFVRegisterAll();
982:   PetscFunctionListFind(PetscFVList, name, &r);
983:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);

985:   if (fvm->ops->destroy) {
986:     (*fvm->ops->destroy)(fvm);
987:     fvm->ops->destroy = NULL;
988:   }
989:   (*r)(fvm);
990:   PetscObjectChangeTypeName((PetscObject) fvm, name);
991:   return(0);
992: }

994: /*@C
995:   PetscFVGetType - Gets the PetscFV type name (as a string) from the object.

997:   Not Collective

999:   Input Parameter:
1000: . fvm  - The PetscFV

1002:   Output Parameter:
1003: . name - The PetscFV type name

1005:   Level: intermediate

1007: .seealso: PetscFVSetType(), PetscFVCreate()
1008: @*/
1009: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1010: {

1016:   PetscFVRegisterAll();
1017:   *name = ((PetscObject) fvm)->type_name;
1018:   return(0);
1019: }

1021: /*@C
1022:   PetscFVView - Views a PetscFV

1024:   Collective on fvm

1026:   Input Parameter:
1027: + fvm - the PetscFV object to view
1028: - v   - the viewer

1030:   Level: beginner

1032: .seealso: PetscFVDestroy()
1033: @*/
1034: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1035: {

1040:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1041:   if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1042:   return(0);
1043: }

1045: /*@
1046:   PetscFVSetFromOptions - sets parameters in a PetscFV from the options database

1048:   Collective on fvm

1050:   Input Parameter:
1051: . fvm - the PetscFV object to set options for

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

1056:   Level: intermediate

1058: .seealso: PetscFVView()
1059: @*/
1060: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1061: {
1062:   const char    *defaultType;
1063:   char           name[256];
1064:   PetscBool      flg;

1069:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1070:   else                                 defaultType = ((PetscObject) fvm)->type_name;
1071:   PetscFVRegisterAll();

1073:   PetscObjectOptionsBegin((PetscObject) fvm);
1074:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1075:   if (flg) {
1076:     PetscFVSetType(fvm, name);
1077:   } else if (!((PetscObject) fvm)->type_name) {
1078:     PetscFVSetType(fvm, defaultType);

1080:   }
1081:   PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1082:   if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1083:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1084:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1085:   PetscLimiterSetFromOptions(fvm->limiter);
1086:   PetscOptionsEnd();
1087:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1088:   return(0);
1089: }

1091: /*@
1092:   PetscFVSetUp - Construct data structures for the PetscFV

1094:   Collective on fvm

1096:   Input Parameter:
1097: . fvm - the PetscFV object to setup

1099:   Level: intermediate

1101: .seealso: PetscFVView(), PetscFVDestroy()
1102: @*/
1103: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1104: {

1109:   PetscLimiterSetUp(fvm->limiter);
1110:   if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1111:   return(0);
1112: }

1114: /*@
1115:   PetscFVDestroy - Destroys a PetscFV object

1117:   Collective on fvm

1119:   Input Parameter:
1120: . fvm - the PetscFV object to destroy

1122:   Level: beginner

1124: .seealso: PetscFVView()
1125: @*/
1126: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1127: {
1128:   PetscInt       i;

1132:   if (!*fvm) return(0);

1135:   if (--((PetscObject)(*fvm))->refct > 0) {*fvm = 0; return(0);}
1136:   ((PetscObject) (*fvm))->refct = 0;

1138:   for (i = 0; i < (*fvm)->numComponents; i++) {
1139:     PetscFree((*fvm)->componentNames[i]);
1140:   }
1141:   PetscFree((*fvm)->componentNames);
1142:   PetscLimiterDestroy(&(*fvm)->limiter);
1143:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1144:   PetscFree((*fvm)->fluxWork);
1145:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1146:   PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);

1148:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1149:   PetscHeaderDestroy(fvm);
1150:   return(0);
1151: }

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

1156:   Collective

1158:   Input Parameter:
1159: . comm - The communicator for the PetscFV object

1161:   Output Parameter:
1162: . fvm - The PetscFV object

1164:   Level: beginner

1166: .seealso: PetscFVSetType(), PETSCFVUPWIND
1167: @*/
1168: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1169: {
1170:   PetscFV        f;

1175:   *fvm = NULL;
1176:   PetscFVInitializePackage();

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

1181:   PetscLimiterCreate(comm, &f->limiter);
1182:   f->numComponents    = 1;
1183:   f->dim              = 0;
1184:   f->computeGradients = PETSC_FALSE;
1185:   f->fluxWork         = NULL;
1186:   PetscCalloc1(f->numComponents,&f->componentNames);

1188:   *fvm = f;
1189:   return(0);
1190: }

1192: /*@
1193:   PetscFVSetLimiter - Set the limiter object

1195:   Logically collective on fvm

1197:   Input Parameters:
1198: + fvm - the PetscFV object
1199: - lim - The PetscLimiter

1201:   Level: intermediate

1203: .seealso: PetscFVGetLimiter()
1204: @*/
1205: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1206: {

1212:   PetscLimiterDestroy(&fvm->limiter);
1213:   PetscObjectReference((PetscObject) lim);
1214:   fvm->limiter = lim;
1215:   return(0);
1216: }

1218: /*@
1219:   PetscFVGetLimiter - Get the limiter object

1221:   Not collective

1223:   Input Parameter:
1224: . fvm - the PetscFV object

1226:   Output Parameter:
1227: . lim - The PetscLimiter

1229:   Level: intermediate

1231: .seealso: PetscFVSetLimiter()
1232: @*/
1233: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1234: {
1238:   *lim = fvm->limiter;
1239:   return(0);
1240: }

1242: /*@
1243:   PetscFVSetNumComponents - Set the number of field components

1245:   Logically collective on fvm

1247:   Input Parameters:
1248: + fvm - the PetscFV object
1249: - comp - The number of components

1251:   Level: intermediate

1253: .seealso: PetscFVGetNumComponents()
1254: @*/
1255: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1256: {

1261:   if (fvm->numComponents != comp) {
1262:     PetscInt i;

1264:     for (i = 0; i < fvm->numComponents; i++) {
1265:       PetscFree(fvm->componentNames[i]);
1266:     }
1267:     PetscFree(fvm->componentNames);
1268:     PetscCalloc1(comp,&fvm->componentNames);
1269:   }
1270:   fvm->numComponents = comp;
1271:   PetscFree(fvm->fluxWork);
1272:   PetscMalloc1(comp, &fvm->fluxWork);
1273:   return(0);
1274: }

1276: /*@
1277:   PetscFVGetNumComponents - Get the number of field components

1279:   Not collective

1281:   Input Parameter:
1282: . fvm - the PetscFV object

1284:   Output Parameter:
1285: , comp - The number of components

1287:   Level: intermediate

1289: .seealso: PetscFVSetNumComponents()
1290: @*/
1291: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1292: {
1296:   *comp = fvm->numComponents;
1297:   return(0);
1298: }

1300: /*@C
1301:   PetscFVSetComponentName - Set the name of a component (used in output and viewing)

1303:   Logically collective on fvm
1304:   Input Parameters:
1305: + fvm - the PetscFV object
1306: . comp - the component number
1307: - name - the component name

1309:   Level: intermediate

1311: .seealso: PetscFVGetComponentName()
1312: @*/
1313: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1314: {

1318:   PetscFree(fvm->componentNames[comp]);
1319:   PetscStrallocpy(name,&fvm->componentNames[comp]);
1320:   return(0);
1321: }

1323: /*@C
1324:   PetscFVGetComponentName - Get the name of a component (used in output and viewing)

1326:   Logically collective on fvm
1327:   Input Parameters:
1328: + fvm - the PetscFV object
1329: - comp - the component number

1331:   Output Parameter:
1332: . name - the component name

1334:   Level: intermediate

1336: .seealso: PetscFVSetComponentName()
1337: @*/
1338: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1339: {
1341:   *name = fvm->componentNames[comp];
1342:   return(0);
1343: }

1345: /*@
1346:   PetscFVSetSpatialDimension - Set the spatial dimension

1348:   Logically collective on fvm

1350:   Input Parameters:
1351: + fvm - the PetscFV object
1352: - dim - The spatial dimension

1354:   Level: intermediate

1356: .seealso: PetscFVGetSpatialDimension()
1357: @*/
1358: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1359: {
1362:   fvm->dim = dim;
1363:   return(0);
1364: }

1366: /*@
1367:   PetscFVGetSpatialDimension - Get the spatial dimension

1369:   Logically collective on fvm

1371:   Input Parameter:
1372: . fvm - the PetscFV object

1374:   Output Parameter:
1375: . dim - The spatial dimension

1377:   Level: intermediate

1379: .seealso: PetscFVSetSpatialDimension()
1380: @*/
1381: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1382: {
1386:   *dim = fvm->dim;
1387:   return(0);
1388: }

1390: /*@
1391:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1393:   Logically collective on fvm

1395:   Input Parameters:
1396: + fvm - the PetscFV object
1397: - computeGradients - Flag to compute cell gradients

1399:   Level: intermediate

1401: .seealso: PetscFVGetComputeGradients()
1402: @*/
1403: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1404: {
1407:   fvm->computeGradients = computeGradients;
1408:   return(0);
1409: }

1411: /*@
1412:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1414:   Not collective

1416:   Input Parameter:
1417: . fvm - the PetscFV object

1419:   Output Parameter:
1420: . computeGradients - Flag to compute cell gradients

1422:   Level: intermediate

1424: .seealso: PetscFVSetComputeGradients()
1425: @*/
1426: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1427: {
1431:   *computeGradients = fvm->computeGradients;
1432:   return(0);
1433: }

1435: /*@
1436:   PetscFVSetQuadrature - Set the quadrature object

1438:   Logically collective on fvm

1440:   Input Parameters:
1441: + fvm - the PetscFV object
1442: - q - The PetscQuadrature

1444:   Level: intermediate

1446: .seealso: PetscFVGetQuadrature()
1447: @*/
1448: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1449: {

1454:   PetscQuadratureDestroy(&fvm->quadrature);
1455:   PetscObjectReference((PetscObject) q);
1456:   fvm->quadrature = q;
1457:   return(0);
1458: }

1460: /*@
1461:   PetscFVGetQuadrature - Get the quadrature object

1463:   Not collective

1465:   Input Parameter:
1466: . fvm - the PetscFV object

1468:   Output Parameter:
1469: . lim - The PetscQuadrature

1471:   Level: intermediate

1473: .seealso: PetscFVSetQuadrature()
1474: @*/
1475: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1476: {
1480:   if (!fvm->quadrature) {
1481:     /* Create default 1-point quadrature */
1482:     PetscReal     *points, *weights;

1485:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1486:     PetscCalloc1(fvm->dim, &points);
1487:     PetscMalloc1(1, &weights);
1488:     weights[0] = 1.0;
1489:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1490:   }
1491:   *q = fvm->quadrature;
1492:   return(0);
1493: }

1495: /*@
1496:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1498:   Not collective

1500:   Input Parameter:
1501: . fvm - The PetscFV object

1503:   Output Parameter:
1504: . sp - The PetscDualSpace object

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

1508:   Level: intermediate

1510: .seealso: PetscFVCreate()
1511: @*/
1512: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1513: {
1517:   if (!fvm->dualSpace) {
1518:     DM              K;
1519:     PetscInt        dim, Nc, c;
1520:     PetscErrorCode  ierr;

1522:     PetscFVGetSpatialDimension(fvm, &dim);
1523:     PetscFVGetNumComponents(fvm, &Nc);
1524:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1525:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1526:     PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1527:     PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1528:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1529:     DMDestroy(&K);
1530:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1531:     /* Should we be using PetscFVGetQuadrature() here? */
1532:     for (c = 0; c < Nc; ++c) {
1533:       PetscQuadrature qc;
1534:       PetscReal      *points, *weights;
1535:       PetscErrorCode  ierr;

1537:       PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1538:       PetscCalloc1(dim, &points);
1539:       PetscCalloc1(Nc, &weights);
1540:       weights[c] = 1.0;
1541:       PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1542:       PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1543:       PetscQuadratureDestroy(&qc);
1544:     }
1545:   }
1546:   *sp = fvm->dualSpace;
1547:   return(0);
1548: }

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

1553:   Not collective

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

1559:   Level: intermediate

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

1563: .seealso: PetscFVCreate()
1564: @*/
1565: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1566: {

1572:   PetscDualSpaceDestroy(&fvm->dualSpace);
1573:   fvm->dualSpace = sp;
1574:   PetscObjectReference((PetscObject) fvm->dualSpace);
1575:   return(0);
1576: }

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

1581:   Not collective

1583:   Input Parameter:
1584: . fvm - The PetscFV object

1586:   Output Parameters:
1587: + B - The basis function values at quadrature points
1588: . D - The basis function derivatives at quadrature points
1589: - H - The basis function second derivatives at quadrature points

1591:   Note:
1592: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1593: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1594: $ 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

1596:   Level: intermediate

1598: .seealso: PetscFEGetDefaultTabulation(), PetscFEGetTabulation(), PetscFERestoreTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1599: @*/
1600: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1601: {
1602:   PetscInt         npoints;
1603:   const PetscReal *points;
1604:   PetscErrorCode   ierr;

1611:   PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1612:   if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1613:   if (B) *B = fvm->B;
1614:   if (D) *D = fvm->D;
1615:   if (H) *H = fvm->H;
1616:   return(0);
1617: }

1619: /*@C
1620:   PetscFVGetTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.

1622:   Not collective

1624:   Input Parameters:
1625: + fvm     - The PetscFV object
1626: . npoints - The number of tabulation points
1627: - points  - The tabulation point coordinates

1629:   Output Parameters:
1630: + B - The basis function values at tabulation points
1631: . D - The basis function derivatives at tabulation points
1632: - H - The basis function second derivatives at tabulation points

1634:   Note:
1635: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1636: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1637: $ 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

1639:   Level: intermediate

1641: .seealso: PetscFEGetTabulation(), PetscFERestoreTabulation(), PetscFEGetDefaultTabulation()
1642: @*/
1643: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1644: {
1645:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1646:   PetscInt         dim;      /* Spatial dimension */
1647:   PetscInt         comp;     /* Field components */
1648:   PetscInt         p, d, c, e;
1649:   PetscErrorCode   ierr;

1657:   PetscFVGetSpatialDimension(fvm, &dim);
1658:   PetscFVGetNumComponents(fvm, &comp);
1659:   if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1660:   if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1661:   if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1662:   if (B) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) (*B)[(p*pdim + d)*comp + c] = 1.0;}
1663:   if (D) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim; ++e) (*D)[((p*pdim + d)*comp + c)*dim + e] = 0.0;}
1664:   if (H) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim*dim; ++e) (*H)[((p*pdim + d)*comp + c)*dim*dim + e] = 0.0;}
1665:   return(0);
1666: }

1668: /*@C
1669:   PetscFVRestoreTabulation - Frees memory from the associated tabulation.

1671:   Not collective

1673:   Input Parameters:
1674: + fvm     - The PetscFV object
1675: . npoints - The number of tabulation points
1676: . points  - The tabulation point coordinates
1677: . B - The basis function values at tabulation points
1678: . D - The basis function derivatives at tabulation points
1679: - H - The basis function second derivatives at tabulation points

1681:   Note:
1682: $ B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1683: $ D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1684: $ 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

1686:   Level: intermediate

1688: .seealso: PetscFVGetTabulation(), PetscFVGetDefaultTabulation()
1689: @*/
1690: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1691: {

1696:   if (B && *B) {PetscFree(*B);}
1697:   if (D && *D) {PetscFree(*D);}
1698:   if (H && *H) {PetscFree(*H);}
1699:   return(0);
1700: }

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

1705:   Input Parameters:
1706: + fvm      - The PetscFV object
1707: . numFaces - The number of cell faces which are not constrained
1708: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1710:   Level: advanced

1712: .seealso: PetscFVCreate()
1713: @*/
1714: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1715: {

1720:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1721:   return(0);
1722: }

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

1727:   Not collective

1729:   Input Parameters:
1730: + fvm          - The PetscFV object for the field being integrated
1731: . prob         - The PetscDS specifing the discretizations and continuum functions
1732: . field        - The field being integrated
1733: . Nf           - The number of faces in the chunk
1734: . fgeom        - The face geometry for each face in the chunk
1735: . neighborVol  - The volume for each pair of cells in the chunk
1736: . uL           - The state from the cell on the left
1737: - uR           - The state from the cell on the right

1739:   Output Parameter
1740: + fluxL        - the left fluxes for each face
1741: - fluxR        - the right fluxes for each face

1743:   Level: developer

1745: .seealso: PetscFVCreate()
1746: @*/
1747: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1748:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1749: {

1754:   if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);}
1755:   return(0);
1756: }

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

1763:   Input Parameter:
1764: . fv - The initial PetscFV

1766:   Output Parameter:
1767: . fvRef - The refined PetscFV

1769:   Level: advanced

1771: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1772: @*/
1773: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1774: {
1775:   PetscDualSpace   Q, Qref;
1776:   DM               K, Kref;
1777:   PetscQuadrature  q, qref;
1778:   CellRefiner      cellRefiner;
1779:   PetscReal       *v0;
1780:   PetscReal       *jac, *invjac;
1781:   PetscInt         numComp, numSubelements, s;
1782:   PetscErrorCode   ierr;

1785:   PetscFVGetDualSpace(fv, &Q);
1786:   PetscFVGetQuadrature(fv, &q);
1787:   PetscDualSpaceGetDM(Q, &K);
1788:   /* Create dual space */
1789:   PetscDualSpaceDuplicate(Q, &Qref);
1790:   DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1791:   PetscDualSpaceSetDM(Qref, Kref);
1792:   DMDestroy(&Kref);
1793:   PetscDualSpaceSetUp(Qref);
1794:   /* Create volume */
1795:   PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1796:   PetscFVSetDualSpace(*fvRef, Qref);
1797:   PetscFVGetNumComponents(fv,    &numComp);
1798:   PetscFVSetNumComponents(*fvRef, numComp);
1799:   PetscFVSetUp(*fvRef);
1800:   /* Create quadrature */
1801:   DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1802:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1803:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1804:   PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1805:   for (s = 0; s < numSubelements; ++s) {
1806:     PetscQuadrature  qs;
1807:     const PetscReal *points, *weights;
1808:     PetscReal       *p, *w;
1809:     PetscInt         dim, Nc, npoints, np;

1811:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1812:     PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1813:     np   = npoints/numSubelements;
1814:     PetscMalloc1(np*dim,&p);
1815:     PetscMalloc1(np*Nc,&w);
1816:     PetscArraycpy(p, &points[s*np*dim], np*dim);
1817:     PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1818:     PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1819:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1820:     PetscQuadratureDestroy(&qs);
1821:   }
1822:   CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1823:   PetscFVSetQuadrature(*fvRef, qref);
1824:   PetscQuadratureDestroy(&qref);
1825:   PetscDualSpaceDestroy(&Qref);
1826:   return(0);
1827: }

1829: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1830: {
1831:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1835:   PetscFree(b);
1836:   return(0);
1837: }

1839: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1840: {
1841:   PetscInt          Nc, c;
1842:   PetscViewerFormat format;
1843:   PetscErrorCode    ierr;

1846:   PetscFVGetNumComponents(fv, &Nc);
1847:   PetscViewerGetFormat(viewer, &format);
1848:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1849:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1850:   for (c = 0; c < Nc; c++) {
1851:     if (fv->componentNames[c]) {
1852:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1853:     }
1854:   }
1855:   return(0);
1856: }

1858: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1859: {
1860:   PetscBool      iascii;

1866:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1867:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1868:   return(0);
1869: }

1871: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1872: {
1874:   return(0);
1875: }

1877: /*
1878:   neighborVol[f*2+0] contains the left  geom
1879:   neighborVol[f*2+1] contains the right geom
1880: */
1881: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1882:                                                          PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1883: {
1884:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1885:   void              *rctx;
1886:   PetscScalar       *flux = fvm->fluxWork;
1887:   const PetscScalar *constants;
1888:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;
1889:   PetscErrorCode     ierr;

1892:   PetscDSGetTotalComponents(prob, &Nc);
1893:   PetscDSGetTotalDimension(prob, &totDim);
1894:   PetscDSGetFieldOffset(prob, field, &off);
1895:   PetscDSGetRiemannSolver(prob, field, &riemann);
1896:   PetscDSGetContext(prob, field, &rctx);
1897:   PetscDSGetConstants(prob, &numConstants, &constants);
1898:   PetscFVGetSpatialDimension(fvm, &dim);
1899:   PetscFVGetNumComponents(fvm, &pdim);
1900:   for (f = 0; f < Nf; ++f) {
1901:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1902:     for (d = 0; d < pdim; ++d) {
1903:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1904:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1905:     }
1906:   }
1907:   return(0);
1908: }

1910: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1911: {
1913:   fvm->ops->setfromoptions          = NULL;
1914:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1915:   fvm->ops->view                    = PetscFVView_Upwind;
1916:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
1917:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
1918:   return(0);
1919: }

1921: /*MC
1922:   PETSCFVUPWIND = "upwind" - A PetscFV object

1924:   Level: intermediate

1926: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1927: M*/

1929: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1930: {
1931:   PetscFV_Upwind *b;
1932:   PetscErrorCode  ierr;

1936:   PetscNewLog(fvm,&b);
1937:   fvm->data = b;

1939:   PetscFVInitialize_Upwind(fvm);
1940:   return(0);
1941: }

1943:  #include <petscblaslapack.h>

1945: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1946: {
1947:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1948:   PetscErrorCode        ierr;

1951:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1952:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1953:   PetscFree(ls);
1954:   return(0);
1955: }

1957: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1958: {
1959:   PetscInt          Nc, c;
1960:   PetscViewerFormat format;
1961:   PetscErrorCode    ierr;

1964:   PetscFVGetNumComponents(fv, &Nc);
1965:   PetscViewerGetFormat(viewer, &format);
1966:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1967:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1968:   for (c = 0; c < Nc; c++) {
1969:     if (fv->componentNames[c]) {
1970:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1971:     }
1972:   }
1973:   return(0);
1974: }

1976: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1977: {
1978:   PetscBool      iascii;

1984:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1985:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1986:   return(0);
1987: }

1989: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1990: {
1992:   return(0);
1993: }

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

2004:   if (debug) {
2005:     PetscMalloc1(m*n,&Aback);
2006:     PetscArraycpy(Aback,A,m*n);
2007:   }

2009:   PetscBLASIntCast(m,&M);
2010:   PetscBLASIntCast(n,&N);
2011:   PetscBLASIntCast(mstride,&lda);
2012:   PetscBLASIntCast(worksize,&ldwork);
2013:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2014:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2015:   PetscFPTrapPop();
2016:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2017:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2019:   /* Extract an explicit representation of Q */
2020:   Q    = Ainv;
2021:   PetscArraycpy(Q,A,mstride*n);
2022:   K    = N;                     /* full rank */
2023:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2024:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

2026:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2027:   Alpha = 1.0;
2028:   ldb   = lda;
2029:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2030:   /* Ainv is Q, overwritten with inverse */

2032:   if (debug) {                      /* Check that pseudo-inverse worked */
2033:     PetscScalar  Beta = 0.0;
2034:     PetscBLASInt ldc;
2035:     K   = N;
2036:     ldc = N;
2037:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2038:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2039:     PetscFree(Aback);
2040:   }
2041:   return(0);
2042: }

2044: /* Overwrites A. Can handle degenerate problems and m<n. */
2045: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2046: {
2047:   PetscBool      debug = PETSC_FALSE;
2048:   PetscScalar   *Brhs, *Aback;
2049:   PetscScalar   *tmpwork;
2050:   PetscReal      rcond;
2051: #if defined (PETSC_USE_COMPLEX)
2052:   PetscInt       rworkSize;
2053:   PetscReal     *rwork;
2054: #endif
2055:   PetscInt       i, j, maxmn;
2056:   PetscBLASInt   M, N, lda, ldb, ldwork;
2057:   PetscBLASInt   nrhs, irank, info;

2061:   if (debug) {
2062:     PetscMalloc1(m*n,&Aback);
2063:     PetscArraycpy(Aback,A,m*n);
2064:   }

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

2074:   PetscBLASIntCast(m,&M);
2075:   PetscBLASIntCast(n,&N);
2076:   PetscBLASIntCast(mstride,&lda);
2077:   PetscBLASIntCast(maxmn,&ldb);
2078:   PetscBLASIntCast(worksize,&ldwork);
2079:   rcond = -1;
2080:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2081:   nrhs  = M;
2082: #if defined(PETSC_USE_COMPLEX)
2083:   rworkSize = 5 * PetscMin(M,N);
2084:   PetscMalloc1(rworkSize,&rwork);
2085:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2086:   PetscFPTrapPop();
2087:   PetscFree(rwork);
2088: #else
2089:   nrhs  = M;
2090:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2091:   PetscFPTrapPop();
2092: #endif
2093:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2094:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2095:   if (irank < PetscMin(M,N)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Rank deficient least squares fit, indicates an isolated cell with two colinear points");

2097:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2098:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2099:   for (i=0; i<n; i++) {
2100:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2101:   }
2102:   return(0);
2103: }

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

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

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

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

2138: /*
2139:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

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

2146:   Level: developer

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

2159:   if (numFaces > maxFaces) {
2160:     if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2161:     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2162:   }
2163:   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) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2169:   else        {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2170:   for (f = 0; f < numFaces; ++f) {
2171:     for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2172:   }
2173:   return(0);
2174: }

2176: /*
2177:   neighborVol[f*2+0] contains the left  geom
2178:   neighborVol[f*2+1] contains the right geom
2179: */
2180: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2181:                                                                PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2182: {
2183:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2184:   void              *rctx;
2185:   PetscScalar       *flux = fvm->fluxWork;
2186:   const PetscScalar *constants;
2187:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;
2188:   PetscErrorCode     ierr;

2191:   PetscDSGetTotalComponents(prob, &Nc);
2192:   PetscDSGetTotalDimension(prob, &totDim);
2193:   PetscDSGetFieldOffset(prob, field, &off);
2194:   PetscDSGetRiemannSolver(prob, field, &riemann);
2195:   PetscDSGetContext(prob, field, &rctx);
2196:   PetscDSGetConstants(prob, &numConstants, &constants);
2197:   PetscFVGetSpatialDimension(fvm, &dim);
2198:   PetscFVGetNumComponents(fvm, &pdim);
2199:   for (f = 0; f < Nf; ++f) {
2200:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2201:     for (d = 0; d < pdim; ++d) {
2202:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2203:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2204:     }
2205:   }
2206:   return(0);
2207: }

2209: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2210: {
2211:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2212:   PetscInt              dim, m, n, nrhs, minwork;
2213:   PetscErrorCode        ierr;

2217:   PetscFVGetSpatialDimension(fvm, &dim);
2218:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2219:   ls->maxFaces = maxFaces;
2220:   m       = ls->maxFaces;
2221:   n       = dim;
2222:   nrhs    = ls->maxFaces;
2223:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2224:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2225:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2226:   return(0);
2227: }

2229: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2230: {
2232:   fvm->ops->setfromoptions          = NULL;
2233:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2234:   fvm->ops->view                    = PetscFVView_LeastSquares;
2235:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2236:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2237:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2238:   return(0);
2239: }

2241: /*MC
2242:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2244:   Level: intermediate

2246: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2247: M*/

2249: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2250: {
2251:   PetscFV_LeastSquares *ls;
2252:   PetscErrorCode        ierr;

2256:   PetscNewLog(fvm, &ls);
2257:   fvm->data = ls;

2259:   ls->maxFaces = -1;
2260:   ls->workSize = -1;
2261:   ls->B        = NULL;
2262:   ls->Binv     = NULL;
2263:   ls->tau      = NULL;
2264:   ls->work     = NULL;

2266:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2267:   PetscFVInitialize_LeastSquares(fvm);
2268:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2269:   return(0);
2270: }

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

2275:   Not collective

2277:   Input parameters:
2278: + fvm      - The PetscFV object
2279: - maxFaces - The maximum number of cell faces

2281:   Level: intermediate

2283: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2284: @*/
2285: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2286: {

2291:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2292:   return(0);
2293: }