Actual source code: fv.c

  1: #include <petsc/private/petscfvimpl.h>
  2: #include <petscdmplex.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:    PetscLimiterViewFromOptions - View from Options

129:    Collective on PetscLimiter

131:    Input Parameters:
132: +  A - the PetscLimiter object to view
133: .  obj - Optional object
134: -  name - command line option

136:    Level: intermediate
137: .seealso:  PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
138: @*/
139: PetscErrorCode  PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
140: {

145:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
146:   return(0);
147: }

149: /*@C
150:   PetscLimiterView - Views a PetscLimiter

152:   Collective on lim

154:   Input Parameter:
155: + lim - the PetscLimiter object to view
156: - v   - the viewer

158:   Level: beginner

160: .seealso: PetscLimiterDestroy()
161: @*/
162: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
163: {

168:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
169:   if (lim->ops->view) {(*lim->ops->view)(lim, v);}
170:   return(0);
171: }

173: /*@
174:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

176:   Collective on lim

178:   Input Parameter:
179: . lim - the PetscLimiter object to set options for

181:   Level: intermediate

183: .seealso: PetscLimiterView()
184: @*/
185: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
186: {
187:   const char    *defaultType;
188:   char           name[256];
189:   PetscBool      flg;

194:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
195:   else                                 defaultType = ((PetscObject) lim)->type_name;
196:   PetscLimiterRegisterAll();

198:   PetscObjectOptionsBegin((PetscObject) lim);
199:   PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
200:   if (flg) {
201:     PetscLimiterSetType(lim, name);
202:   } else if (!((PetscObject) lim)->type_name) {
203:     PetscLimiterSetType(lim, defaultType);
204:   }
205:   if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
206:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
207:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);
208:   PetscOptionsEnd();
209:   PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
210:   return(0);
211: }

213: /*@C
214:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

216:   Collective on lim

218:   Input Parameter:
219: . lim - the PetscLimiter object to setup

221:   Level: intermediate

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

231:   if (lim->ops->setup) {(*lim->ops->setup)(lim);}
232:   return(0);
233: }

235: /*@
236:   PetscLimiterDestroy - Destroys a PetscLimiter object

238:   Collective on lim

240:   Input Parameter:
241: . lim - the PetscLimiter object to destroy

243:   Level: beginner

245: .seealso: PetscLimiterView()
246: @*/
247: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
248: {

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

255:   if (--((PetscObject)(*lim))->refct > 0) {*lim = NULL; return(0);}
256:   ((PetscObject) (*lim))->refct = 0;

258:   if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
259:   PetscHeaderDestroy(lim);
260:   return(0);
261: }

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

266:   Collective

268:   Input Parameter:
269: . comm - The communicator for the PetscLimiter object

271:   Output Parameter:
272: . lim - The PetscLimiter object

274:   Level: beginner

276: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
277: @*/
278: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
279: {
280:   PetscLimiter   l;

285:   PetscCitationsRegister(LimiterCitation,&Limitercite);
286:   *lim = NULL;
287:   PetscFVInitializePackage();

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

291:   *lim = l;
292:   return(0);
293: }

295: /*@
296:   PetscLimiterLimit - Limit the flux

298:   Input Parameters:
299: + lim  - The PetscLimiter
300: - flim - The input field

302:   Output Parameter:
303: . phi  - The limited field

305: Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
306: $ The classical flux-limited formulation is psi(r) where
307: $
308: $ r = (u[0] - u[-1]) / (u[1] - u[0])
309: $
310: $ The second order TVD region is bounded by
311: $
312: $ psi_minmod(r) = min(r,1)      and        psi_superbee(r) = min(2, 2r, max(1,r))
313: $
314: $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
315: $ phi(r)(r+1)/2 in which the reconstructed interface values are
316: $
317: $ u(v) = u[0] + phi(r) (grad u)[0] v
318: $
319: $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
320: $
321: $ 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))
322: $
323: $ For a nicer symmetric formulation, rewrite in terms of
324: $
325: $ f = (u[0] - u[-1]) / (u[1] - u[-1])
326: $
327: $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
328: $
329: $ phi(r) = phi(1/r)
330: $
331: $ becomes
332: $
333: $ w(f) = w(1-f).
334: $
335: $ The limiters below implement this final form w(f). The reference methods are
336: $
337: $ w_minmod(f) = 2 min(f,(1-f))             w_superbee(r) = 4 min((1-f), f)

339:   Level: beginner

341: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
342: @*/
343: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
344: {

350:   (*lim->ops->limit)(lim, flim, phi);
351:   return(0);
352: }

354: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
355: {
356:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
357:   PetscErrorCode    ierr;

360:   PetscFree(l);
361:   return(0);
362: }

364: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
365: {
366:   PetscViewerFormat format;
367:   PetscErrorCode    ierr;

370:   PetscViewerGetFormat(viewer, &format);
371:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
372:   return(0);
373: }

375: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
376: {
377:   PetscBool      iascii;

383:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
384:   if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
385:   return(0);
386: }

388: static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
389: {
391:   *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
392:   return(0);
393: }

395: static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
396: {
398:   lim->ops->view    = PetscLimiterView_Sin;
399:   lim->ops->destroy = PetscLimiterDestroy_Sin;
400:   lim->ops->limit   = PetscLimiterLimit_Sin;
401:   return(0);
402: }

404: /*MC
405:   PETSCLIMITERSIN = "sin" - A PetscLimiter object

407:   Level: intermediate

409: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
410: M*/

412: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
413: {
414:   PetscLimiter_Sin *l;
415:   PetscErrorCode    ierr;

419:   PetscNewLog(lim, &l);
420:   lim->data = l;

422:   PetscLimiterInitialize_Sin(lim);
423:   return(0);
424: }

426: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
427: {
428:   PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
429:   PetscErrorCode    ierr;

432:   PetscFree(l);
433:   return(0);
434: }

436: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
437: {
438:   PetscViewerFormat format;
439:   PetscErrorCode    ierr;

442:   PetscViewerGetFormat(viewer, &format);
443:   PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
444:   return(0);
445: }

447: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
448: {
449:   PetscBool      iascii;

455:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
456:   if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
457:   return(0);
458: }

460: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
461: {
463:   *phi = 0.0;
464:   return(0);
465: }

467: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
468: {
470:   lim->ops->view    = PetscLimiterView_Zero;
471:   lim->ops->destroy = PetscLimiterDestroy_Zero;
472:   lim->ops->limit   = PetscLimiterLimit_Zero;
473:   return(0);
474: }

476: /*MC
477:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

479:   Level: intermediate

481: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
482: M*/

484: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
485: {
486:   PetscLimiter_Zero *l;
487:   PetscErrorCode     ierr;

491:   PetscNewLog(lim, &l);
492:   lim->data = l;

494:   PetscLimiterInitialize_Zero(lim);
495:   return(0);
496: }

498: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
499: {
500:   PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
501:   PetscErrorCode    ierr;

504:   PetscFree(l);
505:   return(0);
506: }

508: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
509: {
510:   PetscViewerFormat format;
511:   PetscErrorCode    ierr;

514:   PetscViewerGetFormat(viewer, &format);
515:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
516:   return(0);
517: }

519: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
520: {
521:   PetscBool      iascii;

527:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
528:   if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
529:   return(0);
530: }

532: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
533: {
535:   *phi = 1.0;
536:   return(0);
537: }

539: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
540: {
542:   lim->ops->view    = PetscLimiterView_None;
543:   lim->ops->destroy = PetscLimiterDestroy_None;
544:   lim->ops->limit   = PetscLimiterLimit_None;
545:   return(0);
546: }

548: /*MC
549:   PETSCLIMITERNONE = "none" - A PetscLimiter object

551:   Level: intermediate

553: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
554: M*/

556: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
557: {
558:   PetscLimiter_None *l;
559:   PetscErrorCode    ierr;

563:   PetscNewLog(lim, &l);
564:   lim->data = l;

566:   PetscLimiterInitialize_None(lim);
567:   return(0);
568: }

570: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
571: {
572:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
573:   PetscErrorCode    ierr;

576:   PetscFree(l);
577:   return(0);
578: }

580: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
581: {
582:   PetscViewerFormat format;
583:   PetscErrorCode    ierr;

586:   PetscViewerGetFormat(viewer, &format);
587:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
588:   return(0);
589: }

591: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
592: {
593:   PetscBool      iascii;

599:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
600:   if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
601:   return(0);
602: }

604: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
605: {
607:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
608:   return(0);
609: }

611: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
612: {
614:   lim->ops->view    = PetscLimiterView_Minmod;
615:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
616:   lim->ops->limit   = PetscLimiterLimit_Minmod;
617:   return(0);
618: }

620: /*MC
621:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

623:   Level: intermediate

625: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
626: M*/

628: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
629: {
630:   PetscLimiter_Minmod *l;
631:   PetscErrorCode    ierr;

635:   PetscNewLog(lim, &l);
636:   lim->data = l;

638:   PetscLimiterInitialize_Minmod(lim);
639:   return(0);
640: }

642: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
643: {
644:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
645:   PetscErrorCode    ierr;

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

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

658:   PetscViewerGetFormat(viewer, &format);
659:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
660:   return(0);
661: }

663: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
664: {
665:   PetscBool      iascii;

671:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
672:   if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
673:   return(0);
674: }

676: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
677: {
679:   *phi = PetscMax(0, 4*f*(1-f));
680:   return(0);
681: }

683: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
684: {
686:   lim->ops->view    = PetscLimiterView_VanLeer;
687:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
688:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
689:   return(0);
690: }

692: /*MC
693:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

695:   Level: intermediate

697: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
698: M*/

700: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
701: {
702:   PetscLimiter_VanLeer *l;
703:   PetscErrorCode    ierr;

707:   PetscNewLog(lim, &l);
708:   lim->data = l;

710:   PetscLimiterInitialize_VanLeer(lim);
711:   return(0);
712: }

714: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
715: {
716:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
717:   PetscErrorCode    ierr;

720:   PetscFree(l);
721:   return(0);
722: }

724: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
725: {
726:   PetscViewerFormat format;
727:   PetscErrorCode    ierr;

730:   PetscViewerGetFormat(viewer, &format);
731:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
732:   return(0);
733: }

735: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
736: {
737:   PetscBool      iascii;

743:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
744:   if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
745:   return(0);
746: }

748: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
749: {
751:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
752:   return(0);
753: }

755: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
756: {
758:   lim->ops->view    = PetscLimiterView_VanAlbada;
759:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
760:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
761:   return(0);
762: }

764: /*MC
765:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

767:   Level: intermediate

769: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
770: M*/

772: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
773: {
774:   PetscLimiter_VanAlbada *l;
775:   PetscErrorCode    ierr;

779:   PetscNewLog(lim, &l);
780:   lim->data = l;

782:   PetscLimiterInitialize_VanAlbada(lim);
783:   return(0);
784: }

786: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
787: {
788:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
789:   PetscErrorCode    ierr;

792:   PetscFree(l);
793:   return(0);
794: }

796: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
797: {
798:   PetscViewerFormat format;
799:   PetscErrorCode    ierr;

802:   PetscViewerGetFormat(viewer, &format);
803:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
804:   return(0);
805: }

807: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
808: {
809:   PetscBool      iascii;

815:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
816:   if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
817:   return(0);
818: }

820: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
821: {
823:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
824:   return(0);
825: }

827: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
828: {
830:   lim->ops->view    = PetscLimiterView_Superbee;
831:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
832:   lim->ops->limit   = PetscLimiterLimit_Superbee;
833:   return(0);
834: }

836: /*MC
837:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

839:   Level: intermediate

841: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
842: M*/

844: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
845: {
846:   PetscLimiter_Superbee *l;
847:   PetscErrorCode    ierr;

851:   PetscNewLog(lim, &l);
852:   lim->data = l;

854:   PetscLimiterInitialize_Superbee(lim);
855:   return(0);
856: }

858: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
859: {
860:   PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
861:   PetscErrorCode    ierr;

864:   PetscFree(l);
865:   return(0);
866: }

868: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
869: {
870:   PetscViewerFormat format;
871:   PetscErrorCode    ierr;

874:   PetscViewerGetFormat(viewer, &format);
875:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
876:   return(0);
877: }

879: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
880: {
881:   PetscBool      iascii;

887:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
888:   if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
889:   return(0);
890: }

892: /* aka Barth-Jespersen */
893: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
894: {
896:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
897:   return(0);
898: }

900: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
901: {
903:   lim->ops->view    = PetscLimiterView_MC;
904:   lim->ops->destroy = PetscLimiterDestroy_MC;
905:   lim->ops->limit   = PetscLimiterLimit_MC;
906:   return(0);
907: }

909: /*MC
910:   PETSCLIMITERMC = "mc" - A PetscLimiter object

912:   Level: intermediate

914: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
915: M*/

917: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
918: {
919:   PetscLimiter_MC *l;
920:   PetscErrorCode    ierr;

924:   PetscNewLog(lim, &l);
925:   lim->data = l;

927:   PetscLimiterInitialize_MC(lim);
928:   return(0);
929: }

931: PetscClassId PETSCFV_CLASSID = 0;

933: PetscFunctionList PetscFVList              = NULL;
934: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

936: /*@C
937:   PetscFVRegister - Adds a new PetscFV implementation

939:   Not Collective

941:   Input Parameters:
942: + name        - The name of a new user-defined creation routine
943: - create_func - The creation routine itself

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

948:   Sample usage:
949: .vb
950:     PetscFVRegister("my_fv", MyPetscFVCreate);
951: .ve

953:   Then, your PetscFV type can be chosen with the procedural interface via
954: .vb
955:     PetscFVCreate(MPI_Comm, PetscFV *);
956:     PetscFVSetType(PetscFV, "my_fv");
957: .ve
958:    or at runtime via the option
959: .vb
960:     -petscfv_type my_fv
961: .ve

963:   Level: advanced

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

967: @*/
968: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
969: {

973:   PetscFunctionListAdd(&PetscFVList, sname, function);
974:   return(0);
975: }

977: /*@C
978:   PetscFVSetType - Builds a particular PetscFV

980:   Collective on fvm

982:   Input Parameters:
983: + fvm  - The PetscFV object
984: - name - The kind of FVM space

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

989:   Level: intermediate

991: .seealso: PetscFVGetType(), PetscFVCreate()
992: @*/
993: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
994: {
995:   PetscErrorCode (*r)(PetscFV);
996:   PetscBool      match;

1001:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1002:   if (match) return(0);

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

1008:   if (fvm->ops->destroy) {
1009:     (*fvm->ops->destroy)(fvm);
1010:     fvm->ops->destroy = NULL;
1011:   }
1012:   (*r)(fvm);
1013:   PetscObjectChangeTypeName((PetscObject) fvm, name);
1014:   return(0);
1015: }

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

1020:   Not Collective

1022:   Input Parameter:
1023: . fvm  - The PetscFV

1025:   Output Parameter:
1026: . name - The PetscFV type name

1028:   Level: intermediate

1030: .seealso: PetscFVSetType(), PetscFVCreate()
1031: @*/
1032: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1033: {

1039:   PetscFVRegisterAll();
1040:   *name = ((PetscObject) fvm)->type_name;
1041:   return(0);
1042: }

1044: /*@C
1045:    PetscFVViewFromOptions - View from Options

1047:    Collective on PetscFV

1049:    Input Parameters:
1050: +  A - the PetscFV object
1051: .  obj - Optional object
1052: -  name - command line option

1054:    Level: intermediate
1055: .seealso:  PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
1056: @*/
1057: PetscErrorCode  PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
1058: {

1063:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
1064:   return(0);
1065: }

1067: /*@C
1068:   PetscFVView - Views a PetscFV

1070:   Collective on fvm

1072:   Input Parameter:
1073: + fvm - the PetscFV object to view
1074: - v   - the viewer

1076:   Level: beginner

1078: .seealso: PetscFVDestroy()
1079: @*/
1080: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1081: {

1086:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1087:   if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1088:   return(0);
1089: }

1091: /*@
1092:   PetscFVSetFromOptions - sets parameters in a PetscFV from the options database

1094:   Collective on fvm

1096:   Input Parameter:
1097: . fvm - the PetscFV object to set options for

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

1102:   Level: intermediate

1104: .seealso: PetscFVView()
1105: @*/
1106: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1107: {
1108:   const char    *defaultType;
1109:   char           name[256];
1110:   PetscBool      flg;

1115:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1116:   else                                 defaultType = ((PetscObject) fvm)->type_name;
1117:   PetscFVRegisterAll();

1119:   PetscObjectOptionsBegin((PetscObject) fvm);
1120:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1121:   if (flg) {
1122:     PetscFVSetType(fvm, name);
1123:   } else if (!((PetscObject) fvm)->type_name) {
1124:     PetscFVSetType(fvm, defaultType);

1126:   }
1127:   PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1128:   if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1129:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1130:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1131:   PetscLimiterSetFromOptions(fvm->limiter);
1132:   PetscOptionsEnd();
1133:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1134:   return(0);
1135: }

1137: /*@
1138:   PetscFVSetUp - Construct data structures for the PetscFV

1140:   Collective on fvm

1142:   Input Parameter:
1143: . fvm - the PetscFV object to setup

1145:   Level: intermediate

1147: .seealso: PetscFVView(), PetscFVDestroy()
1148: @*/
1149: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1150: {

1155:   PetscLimiterSetUp(fvm->limiter);
1156:   if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1157:   return(0);
1158: }

1160: /*@
1161:   PetscFVDestroy - Destroys a PetscFV object

1163:   Collective on fvm

1165:   Input Parameter:
1166: . fvm - the PetscFV object to destroy

1168:   Level: beginner

1170: .seealso: PetscFVView()
1171: @*/
1172: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1173: {
1174:   PetscInt       i;

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

1181:   if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; return(0);}
1182:   ((PetscObject) (*fvm))->refct = 0;

1184:   for (i = 0; i < (*fvm)->numComponents; i++) {
1185:     PetscFree((*fvm)->componentNames[i]);
1186:   }
1187:   PetscFree((*fvm)->componentNames);
1188:   PetscLimiterDestroy(&(*fvm)->limiter);
1189:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1190:   PetscFree((*fvm)->fluxWork);
1191:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1192:   PetscTabulationDestroy(&(*fvm)->T);

1194:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1195:   PetscHeaderDestroy(fvm);
1196:   return(0);
1197: }

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

1202:   Collective

1204:   Input Parameter:
1205: . comm - The communicator for the PetscFV object

1207:   Output Parameter:
1208: . fvm - The PetscFV object

1210:   Level: beginner

1212: .seealso: PetscFVSetType(), PETSCFVUPWIND
1213: @*/
1214: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1215: {
1216:   PetscFV        f;

1221:   *fvm = NULL;
1222:   PetscFVInitializePackage();

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

1227:   PetscLimiterCreate(comm, &f->limiter);
1228:   f->numComponents    = 1;
1229:   f->dim              = 0;
1230:   f->computeGradients = PETSC_FALSE;
1231:   f->fluxWork         = NULL;
1232:   PetscCalloc1(f->numComponents,&f->componentNames);

1234:   *fvm = f;
1235:   return(0);
1236: }

1238: /*@
1239:   PetscFVSetLimiter - Set the limiter object

1241:   Logically collective on fvm

1243:   Input Parameters:
1244: + fvm - the PetscFV object
1245: - lim - The PetscLimiter

1247:   Level: intermediate

1249: .seealso: PetscFVGetLimiter()
1250: @*/
1251: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1252: {

1258:   PetscLimiterDestroy(&fvm->limiter);
1259:   PetscObjectReference((PetscObject) lim);
1260:   fvm->limiter = lim;
1261:   return(0);
1262: }

1264: /*@
1265:   PetscFVGetLimiter - Get the limiter object

1267:   Not collective

1269:   Input Parameter:
1270: . fvm - the PetscFV object

1272:   Output Parameter:
1273: . lim - The PetscLimiter

1275:   Level: intermediate

1277: .seealso: PetscFVSetLimiter()
1278: @*/
1279: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1280: {
1284:   *lim = fvm->limiter;
1285:   return(0);
1286: }

1288: /*@
1289:   PetscFVSetNumComponents - Set the number of field components

1291:   Logically collective on fvm

1293:   Input Parameters:
1294: + fvm - the PetscFV object
1295: - comp - The number of components

1297:   Level: intermediate

1299: .seealso: PetscFVGetNumComponents()
1300: @*/
1301: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1302: {

1307:   if (fvm->numComponents != comp) {
1308:     PetscInt i;

1310:     for (i = 0; i < fvm->numComponents; i++) {
1311:       PetscFree(fvm->componentNames[i]);
1312:     }
1313:     PetscFree(fvm->componentNames);
1314:     PetscCalloc1(comp,&fvm->componentNames);
1315:   }
1316:   fvm->numComponents = comp;
1317:   PetscFree(fvm->fluxWork);
1318:   PetscMalloc1(comp, &fvm->fluxWork);
1319:   return(0);
1320: }

1322: /*@
1323:   PetscFVGetNumComponents - Get the number of field components

1325:   Not collective

1327:   Input Parameter:
1328: . fvm - the PetscFV object

1330:   Output Parameter:
1331: , comp - The number of components

1333:   Level: intermediate

1335: .seealso: PetscFVSetNumComponents()
1336: @*/
1337: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1338: {
1342:   *comp = fvm->numComponents;
1343:   return(0);
1344: }

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

1349:   Logically collective on fvm
1350:   Input Parameters:
1351: + fvm - the PetscFV object
1352: . comp - the component number
1353: - name - the component name

1355:   Level: intermediate

1357: .seealso: PetscFVGetComponentName()
1358: @*/
1359: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1360: {

1364:   PetscFree(fvm->componentNames[comp]);
1365:   PetscStrallocpy(name,&fvm->componentNames[comp]);
1366:   return(0);
1367: }

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

1372:   Logically collective on fvm
1373:   Input Parameters:
1374: + fvm - the PetscFV object
1375: - comp - the component number

1377:   Output Parameter:
1378: . name - the component name

1380:   Level: intermediate

1382: .seealso: PetscFVSetComponentName()
1383: @*/
1384: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1385: {
1387:   *name = fvm->componentNames[comp];
1388:   return(0);
1389: }

1391: /*@
1392:   PetscFVSetSpatialDimension - Set the spatial dimension

1394:   Logically collective on fvm

1396:   Input Parameters:
1397: + fvm - the PetscFV object
1398: - dim - The spatial dimension

1400:   Level: intermediate

1402: .seealso: PetscFVGetSpatialDimension()
1403: @*/
1404: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1405: {
1408:   fvm->dim = dim;
1409:   return(0);
1410: }

1412: /*@
1413:   PetscFVGetSpatialDimension - Get the spatial dimension

1415:   Logically collective on fvm

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

1420:   Output Parameter:
1421: . dim - The spatial dimension

1423:   Level: intermediate

1425: .seealso: PetscFVSetSpatialDimension()
1426: @*/
1427: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1428: {
1432:   *dim = fvm->dim;
1433:   return(0);
1434: }

1436: /*@
1437:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1439:   Logically collective on fvm

1441:   Input Parameters:
1442: + fvm - the PetscFV object
1443: - computeGradients - Flag to compute cell gradients

1445:   Level: intermediate

1447: .seealso: PetscFVGetComputeGradients()
1448: @*/
1449: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1450: {
1453:   fvm->computeGradients = computeGradients;
1454:   return(0);
1455: }

1457: /*@
1458:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1460:   Not collective

1462:   Input Parameter:
1463: . fvm - the PetscFV object

1465:   Output Parameter:
1466: . computeGradients - Flag to compute cell gradients

1468:   Level: intermediate

1470: .seealso: PetscFVSetComputeGradients()
1471: @*/
1472: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1473: {
1477:   *computeGradients = fvm->computeGradients;
1478:   return(0);
1479: }

1481: /*@
1482:   PetscFVSetQuadrature - Set the quadrature object

1484:   Logically collective on fvm

1486:   Input Parameters:
1487: + fvm - the PetscFV object
1488: - q - The PetscQuadrature

1490:   Level: intermediate

1492: .seealso: PetscFVGetQuadrature()
1493: @*/
1494: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1495: {

1500:   PetscQuadratureDestroy(&fvm->quadrature);
1501:   PetscObjectReference((PetscObject) q);
1502:   fvm->quadrature = q;
1503:   return(0);
1504: }

1506: /*@
1507:   PetscFVGetQuadrature - Get the quadrature object

1509:   Not collective

1511:   Input Parameter:
1512: . fvm - the PetscFV object

1514:   Output Parameter:
1515: . lim - The PetscQuadrature

1517:   Level: intermediate

1519: .seealso: PetscFVSetQuadrature()
1520: @*/
1521: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1522: {
1526:   if (!fvm->quadrature) {
1527:     /* Create default 1-point quadrature */
1528:     PetscReal     *points, *weights;

1531:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1532:     PetscCalloc1(fvm->dim, &points);
1533:     PetscMalloc1(1, &weights);
1534:     weights[0] = 1.0;
1535:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1536:   }
1537:   *q = fvm->quadrature;
1538:   return(0);
1539: }

1541: /*@
1542:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1544:   Not collective

1546:   Input Parameter:
1547: . fvm - The PetscFV object

1549:   Output Parameter:
1550: . sp - The PetscDualSpace object

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

1554:   Level: intermediate

1556: .seealso: PetscFVCreate()
1557: @*/
1558: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1559: {
1563:   if (!fvm->dualSpace) {
1564:     DM              K;
1565:     PetscInt        dim, Nc, c;
1566:     PetscErrorCode  ierr;

1568:     PetscFVGetSpatialDimension(fvm, &dim);
1569:     PetscFVGetNumComponents(fvm, &Nc);
1570:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1571:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1572:     PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1573:     PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1574:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1575:     DMDestroy(&K);
1576:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1577:     /* Should we be using PetscFVGetQuadrature() here? */
1578:     for (c = 0; c < Nc; ++c) {
1579:       PetscQuadrature qc;
1580:       PetscReal      *points, *weights;
1581:       PetscErrorCode  ierr;

1583:       PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1584:       PetscCalloc1(dim, &points);
1585:       PetscCalloc1(Nc, &weights);
1586:       weights[c] = 1.0;
1587:       PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1588:       PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1589:       PetscQuadratureDestroy(&qc);
1590:     }
1591:     PetscDualSpaceSetUp(fvm->dualSpace);
1592:   }
1593:   *sp = fvm->dualSpace;
1594:   return(0);
1595: }

1597: /*@
1598:   PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product

1600:   Not collective

1602:   Input Parameters:
1603: + fvm - The PetscFV object
1604: - sp  - The PetscDualSpace object

1606:   Level: intermediate

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

1610: .seealso: PetscFVCreate()
1611: @*/
1612: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1613: {

1619:   PetscDualSpaceDestroy(&fvm->dualSpace);
1620:   fvm->dualSpace = sp;
1621:   PetscObjectReference((PetscObject) fvm->dualSpace);
1622:   return(0);
1623: }

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

1628:   Not collective

1630:   Input Parameter:
1631: . fvm - The PetscFV object

1633:   Output Parameter:
1634: . T - The basis function values and derviatives at quadrature points

1636:   Note:
1637: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1638: $ 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
1639: $ 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

1641:   Level: intermediate

1643: .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1644: @*/
1645: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1646: {
1647:   PetscInt         npoints;
1648:   const PetscReal *points;
1649:   PetscErrorCode   ierr;

1654:   PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1655:   if (!fvm->T) {PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);}
1656:   *T = fvm->T;
1657:   return(0);
1658: }

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

1663:   Not collective

1665:   Input Parameters:
1666: + fvm     - The PetscFV object
1667: . nrepl   - The number of replicas
1668: . npoints - The number of tabulation points in a replica
1669: . points  - The tabulation point coordinates
1670: - K       - The order of derivative to tabulate

1672:   Output Parameter:
1673: . T - The basis function values and derviative at tabulation points

1675:   Note:
1676: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1677: $ 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
1678: $ 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

1680:   Level: intermediate

1682: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1683: @*/
1684: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1685: {
1686:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1687:   PetscInt         cdim;     /* Spatial dimension */
1688:   PetscInt         Nc;       /* Field components */
1689:   PetscInt         k, p, d, c, e;
1690:   PetscErrorCode   ierr;

1693:   if (!npoints || K < 0) {
1694:     *T = NULL;
1695:     return(0);
1696:   }
1700:   PetscFVGetSpatialDimension(fvm, &cdim);
1701:   PetscFVGetNumComponents(fvm, &Nc);
1702:   PetscMalloc1(1, T);
1703:   (*T)->K    = !cdim ? 0 : K;
1704:   (*T)->Nr   = nrepl;
1705:   (*T)->Np   = npoints;
1706:   (*T)->Nb   = pdim;
1707:   (*T)->Nc   = Nc;
1708:   (*T)->cdim = cdim;
1709:   PetscMalloc1((*T)->K+1, &(*T)->T);
1710:   for (k = 0; k <= (*T)->K; ++k) {
1711:     PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
1712:   }
1713:   if (K >= 0) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) (*T)->T[0][(p*pdim + d)*Nc + c] = 1.0;}
1714:   if (K >= 1) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim; ++e) (*T)->T[1][((p*pdim + d)*Nc + c)*cdim + e] = 0.0;}
1715:   if (K >= 2) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim*cdim; ++e) (*T)->T[2][((p*pdim + d)*Nc + c)*cdim*cdim + e] = 0.0;}
1716:   return(0);
1717: }

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

1722:   Input Parameters:
1723: + fvm      - The PetscFV object
1724: . numFaces - The number of cell faces which are not constrained
1725: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1727:   Level: advanced

1729: .seealso: PetscFVCreate()
1730: @*/
1731: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1732: {

1737:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1738:   return(0);
1739: }

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

1744:   Not collective

1746:   Input Parameters:
1747: + fvm          - The PetscFV object for the field being integrated
1748: . prob         - The PetscDS specifing the discretizations and continuum functions
1749: . field        - The field being integrated
1750: . Nf           - The number of faces in the chunk
1751: . fgeom        - The face geometry for each face in the chunk
1752: . neighborVol  - The volume for each pair of cells in the chunk
1753: . uL           - The state from the cell on the left
1754: - uR           - The state from the cell on the right

1756:   Output Parameter:
1757: + fluxL        - the left fluxes for each face
1758: - fluxR        - the right fluxes for each face

1760:   Level: developer

1762: .seealso: PetscFVCreate()
1763: @*/
1764: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1765:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1766: {

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

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

1780:   Input Parameter:
1781: . fv - The initial PetscFV

1783:   Output Parameter:
1784: . fvRef - The refined PetscFV

1786:   Level: advanced

1788: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1789: @*/
1790: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1791: {
1792:   PetscDualSpace    Q, Qref;
1793:   DM                K, Kref;
1794:   PetscQuadrature   q, qref;
1795:   DMPolytopeType    ct;
1796:   DMPlexCellRefiner cr;
1797:   PetscReal        *v0;
1798:   PetscReal        *jac, *invjac;
1799:   PetscInt          numComp, numSubelements, s;
1800:   PetscErrorCode    ierr;

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

1830:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1831:     PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1832:     np   = npoints/numSubelements;
1833:     PetscMalloc1(np*dim,&p);
1834:     PetscMalloc1(np*Nc,&w);
1835:     PetscArraycpy(p, &points[s*np*dim], np*dim);
1836:     PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1837:     PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1838:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1839:     PetscQuadratureDestroy(&qs);
1840:   }
1841:   PetscFVSetQuadrature(*fvRef, qref);
1842:   DMPlexCellRefinerDestroy(&cr);
1843:   PetscQuadratureDestroy(&qref);
1844:   PetscDualSpaceDestroy(&Qref);
1845:   return(0);
1846: }

1848: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1849: {
1850:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1854:   PetscFree(b);
1855:   return(0);
1856: }

1858: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1859: {
1860:   PetscInt          Nc, c;
1861:   PetscViewerFormat format;
1862:   PetscErrorCode    ierr;

1865:   PetscFVGetNumComponents(fv, &Nc);
1866:   PetscViewerGetFormat(viewer, &format);
1867:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1868:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1869:   for (c = 0; c < Nc; c++) {
1870:     if (fv->componentNames[c]) {
1871:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1872:     }
1873:   }
1874:   return(0);
1875: }

1877: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1878: {
1879:   PetscBool      iascii;

1885:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1886:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1887:   return(0);
1888: }

1890: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1891: {
1893:   return(0);
1894: }

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

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

1929: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1930: {
1932:   fvm->ops->setfromoptions          = NULL;
1933:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1934:   fvm->ops->view                    = PetscFVView_Upwind;
1935:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
1936:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
1937:   return(0);
1938: }

1940: /*MC
1941:   PETSCFVUPWIND = "upwind" - A PetscFV object

1943:   Level: intermediate

1945: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1946: M*/

1948: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1949: {
1950:   PetscFV_Upwind *b;
1951:   PetscErrorCode  ierr;

1955:   PetscNewLog(fvm,&b);
1956:   fvm->data = b;

1958:   PetscFVInitialize_Upwind(fvm);
1959:   return(0);
1960: }

1962: #include <petscblaslapack.h>

1964: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1965: {
1966:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1967:   PetscErrorCode        ierr;

1970:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1971:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1972:   PetscFree(ls);
1973:   return(0);
1974: }

1976: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1977: {
1978:   PetscInt          Nc, c;
1979:   PetscViewerFormat format;
1980:   PetscErrorCode    ierr;

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

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

2003:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2004:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2005:   return(0);
2006: }

2008: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2009: {
2011:   return(0);
2012: }

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

2023:   if (debug) {
2024:     PetscMalloc1(m*n,&Aback);
2025:     PetscArraycpy(Aback,A,m*n);
2026:   }

2028:   PetscBLASIntCast(m,&M);
2029:   PetscBLASIntCast(n,&N);
2030:   PetscBLASIntCast(mstride,&lda);
2031:   PetscBLASIntCast(worksize,&ldwork);
2032:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2033:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
2034:   PetscFPTrapPop();
2035:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2036:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2038:   /* Extract an explicit representation of Q */
2039:   Q    = Ainv;
2040:   PetscArraycpy(Q,A,mstride*n);
2041:   K    = N;                     /* full rank */
2042:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2043:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

2045:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2046:   Alpha = 1.0;
2047:   ldb   = lda;
2048:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2049:   /* Ainv is Q, overwritten with inverse */

2051:   if (debug) {                      /* Check that pseudo-inverse worked */
2052:     PetscScalar  Beta = 0.0;
2053:     PetscBLASInt ldc;
2054:     K   = N;
2055:     ldc = N;
2056:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2057:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2058:     PetscFree(Aback);
2059:   }
2060:   return(0);
2061: }

2063: /* Overwrites A. Can handle degenerate problems and m<n. */
2064: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2065: {
2066:   PetscBool      debug = PETSC_FALSE;
2067:   PetscScalar   *Brhs, *Aback;
2068:   PetscScalar   *tmpwork;
2069:   PetscReal      rcond;
2070: #if defined (PETSC_USE_COMPLEX)
2071:   PetscInt       rworkSize;
2072:   PetscReal     *rwork;
2073: #endif
2074:   PetscInt       i, j, maxmn;
2075:   PetscBLASInt   M, N, lda, ldb, ldwork;
2076:   PetscBLASInt   nrhs, irank, info;

2080:   if (debug) {
2081:     PetscMalloc1(m*n,&Aback);
2082:     PetscArraycpy(Aback,A,m*n);
2083:   }

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

2093:   PetscBLASIntCast(m,&M);
2094:   PetscBLASIntCast(n,&N);
2095:   PetscBLASIntCast(mstride,&lda);
2096:   PetscBLASIntCast(maxmn,&ldb);
2097:   PetscBLASIntCast(worksize,&ldwork);
2098:   rcond = -1;
2099:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2100:   nrhs  = M;
2101: #if defined(PETSC_USE_COMPLEX)
2102:   rworkSize = 5 * PetscMin(M,N);
2103:   PetscMalloc1(rworkSize,&rwork);
2104:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info));
2105:   PetscFPTrapPop();
2106:   PetscFree(rwork);
2107: #else
2108:   nrhs  = M;
2109:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info));
2110:   PetscFPTrapPop();
2111: #endif
2112:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2113:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2114:   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");
2115:   return(0);
2116: }

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

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

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

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

2151: /*
2152:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

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

2159:   Level: developer

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

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

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

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

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

2239:   PetscFVGetSpatialDimension(fvm, &dim);
2240:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2241:   ls->maxFaces = maxFaces;
2242:   m       = ls->maxFaces;
2243:   n       = dim;
2244:   nrhs    = ls->maxFaces;
2245:   minmn   = PetscMin(m,n);
2246:   maxmn   = PetscMax(m,n);
2247:   ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2248:   PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);
2249:   return(0);
2250: }

2252: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2253: {
2255:   fvm->ops->setfromoptions          = NULL;
2256:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2257:   fvm->ops->view                    = PetscFVView_LeastSquares;
2258:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2259:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2260:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2261:   return(0);
2262: }

2264: /*MC
2265:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2267:   Level: intermediate

2269: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2270: M*/

2272: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2273: {
2274:   PetscFV_LeastSquares *ls;
2275:   PetscErrorCode        ierr;

2279:   PetscNewLog(fvm, &ls);
2280:   fvm->data = ls;

2282:   ls->maxFaces = -1;
2283:   ls->workSize = -1;
2284:   ls->B        = NULL;
2285:   ls->Binv     = NULL;
2286:   ls->tau      = NULL;
2287:   ls->work     = NULL;

2289:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2290:   PetscFVInitialize_LeastSquares(fvm);
2291:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2292:   return(0);
2293: }

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

2298:   Not collective

2300:   Input parameters:
2301: + fvm      - The PetscFV object
2302: - maxFaces - The maximum number of cell faces

2304:   Level: intermediate

2306: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2307: @*/
2308: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2309: {

2314:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2315:   return(0);
2316: }