Actual source code: fv.c

petsc-master 2020-01-22
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:    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 = 0; 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 = 0; 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:   }
1592:   *sp = fvm->dualSpace;
1593:   return(0);
1594: }

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

1599:   Not collective

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

1605:   Level: intermediate

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

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

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

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

1627:   Not collective

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

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

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

1640:   Level: intermediate

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

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

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

1662:   Not collective

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

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

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

1679:   Level: intermediate

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

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

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

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

1726:   Level: advanced

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

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

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

1743:   Not collective

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

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

1759:   Level: developer

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

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

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

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

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

1785:   Level: advanced

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

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

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

1845: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1846: {
1847:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1851:   PetscFree(b);
1852:   return(0);
1853: }

1855: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1856: {
1857:   PetscInt          Nc, c;
1858:   PetscViewerFormat format;
1859:   PetscErrorCode    ierr;

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

1874: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1875: {
1876:   PetscBool      iascii;

1882:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1883:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1884:   return(0);
1885: }

1887: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1888: {
1890:   return(0);
1891: }

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

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

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

1937: /*MC
1938:   PETSCFVUPWIND = "upwind" - A PetscFV object

1940:   Level: intermediate

1942: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1943: M*/

1945: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1946: {
1947:   PetscFV_Upwind *b;
1948:   PetscErrorCode  ierr;

1952:   PetscNewLog(fvm,&b);
1953:   fvm->data = b;

1955:   PetscFVInitialize_Upwind(fvm);
1956:   return(0);
1957: }

1959:  #include <petscblaslapack.h>

1961: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1962: {
1963:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1964:   PetscErrorCode        ierr;

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

1973: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1974: {
1975:   PetscInt          Nc, c;
1976:   PetscViewerFormat format;
1977:   PetscErrorCode    ierr;

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

1992: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1993: {
1994:   PetscBool      iascii;

2000:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2001:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2002:   return(0);
2003: }

2005: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2006: {
2008:   return(0);
2009: }

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

2020:   if (debug) {
2021:     PetscMalloc1(m*n,&Aback);
2022:     PetscArraycpy(Aback,A,m*n);
2023:   }

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

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

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

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

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

2077:   if (debug) {
2078:     PetscMalloc1(m*n,&Aback);
2079:     PetscArraycpy(Aback,A,m*n);
2080:   }

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

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

2113:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2114:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2115:   for (i=0; i<n; i++) {
2116:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2117:   }
2118:   return(0);
2119: }

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

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

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

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

2154: /*
2155:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

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

2162:   Level: developer

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

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

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

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

2225: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2226: {
2227:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2228:   PetscInt              dim, m, n, nrhs, minwork;
2229:   PetscErrorCode        ierr;

2233:   PetscFVGetSpatialDimension(fvm, &dim);
2234:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2235:   ls->maxFaces = maxFaces;
2236:   m       = ls->maxFaces;
2237:   n       = dim;
2238:   nrhs    = ls->maxFaces;
2239:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2240:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2241:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2242:   return(0);
2243: }

2245: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2246: {
2248:   fvm->ops->setfromoptions          = NULL;
2249:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2250:   fvm->ops->view                    = PetscFVView_LeastSquares;
2251:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2252:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2253:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2254:   return(0);
2255: }

2257: /*MC
2258:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2260:   Level: intermediate

2262: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2263: M*/

2265: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2266: {
2267:   PetscFV_LeastSquares *ls;
2268:   PetscErrorCode        ierr;

2272:   PetscNewLog(fvm, &ls);
2273:   fvm->data = ls;

2275:   ls->maxFaces = -1;
2276:   ls->workSize = -1;
2277:   ls->B        = NULL;
2278:   ls->Binv     = NULL;
2279:   ls->tau      = NULL;
2280:   ls->work     = NULL;

2282:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2283:   PetscFVInitialize_LeastSquares(fvm);
2284:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2285:   return(0);
2286: }

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

2291:   Not collective

2293:   Input parameters:
2294: + fvm      - The PetscFV object
2295: - maxFaces - The maximum number of cell faces

2297:   Level: intermediate

2299: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2300: @*/
2301: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2302: {

2307:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2308:   return(0);
2309: }