Actual source code: dtfv.c

petsc-3.9.4 2018-09-11
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: .keywords: PetscLimiter, register
 48: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()

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

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

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

 63:   Collective on PetscLimiter

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

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

 72:   Level: intermediate

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

 85:   PetscObjectTypeCompare((PetscObject) lim, name, &match);
 86:   if (match) return(0);

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

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

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

104:   Not Collective

106:   Input Parameter:
107: . lim  - The PetscLimiter

109:   Output Parameter:
110: . name - The PetscLimiter type name

112:   Level: intermediate

114: .keywords: PetscLimiter, get, type, name
115: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
116: @*/
117: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
118: {

124:   PetscLimiterRegisterAll();
125:   *name = ((PetscObject) lim)->type_name;
126:   return(0);
127: }

129: /*@C
130:   PetscLimiterView - Views a PetscLimiter

132:   Collective on PetscLimiter

134:   Input Parameter:
135: + lim - the PetscLimiter object to view
136: - v   - the viewer

138:   Level: developer

140: .seealso: PetscLimiterDestroy()
141: @*/
142: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
143: {

148:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
149:   if (lim->ops->view) {(*lim->ops->view)(lim, v);}
150:   return(0);
151: }

153: /*@
154:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

156:   Collective on PetscLimiter

158:   Input Parameter:
159: . lim - the PetscLimiter object to set options for

161:   Level: developer

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

174:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
175:   else                                 defaultType = ((PetscObject) lim)->type_name;
176:   PetscLimiterRegisterAll();

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

193: /*@C
194:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

196:   Collective on PetscLimiter

198:   Input Parameter:
199: . lim - the PetscLimiter object to setup

201:   Level: developer

203: .seealso: PetscLimiterView(), PetscLimiterDestroy()
204: @*/
205: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
206: {

211:   if (lim->ops->setup) {(*lim->ops->setup)(lim);}
212:   return(0);
213: }

215: /*@
216:   PetscLimiterDestroy - Destroys a PetscLimiter object

218:   Collective on PetscLimiter

220:   Input Parameter:
221: . lim - the PetscLimiter object to destroy

223:   Level: developer

225: .seealso: PetscLimiterView()
226: @*/
227: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
228: {

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

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

238:   if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
239:   PetscHeaderDestroy(lim);
240:   return(0);
241: }

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

246:   Collective on MPI_Comm

248:   Input Parameter:
249: . comm - The communicator for the PetscLimiter object

251:   Output Parameter:
252: . lim - The PetscLimiter object

254:   Level: beginner

256: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
257: @*/
258: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
259: {
260:   PetscLimiter   l;

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

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

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

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

317:   (*lim->ops->limit)(lim, flim, phi);
318:   return(0);
319: }

321: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
322: {
323:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
324:   PetscErrorCode    ierr;

327:   PetscFree(l);
328:   return(0);
329: }

331: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
332: {
333:   PetscViewerFormat format;
334:   PetscErrorCode    ierr;

337:   PetscViewerGetFormat(viewer, &format);
338:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
339:   return(0);
340: }

342: PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
343: {
344:   PetscBool      iascii;

350:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
351:   if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
352:   return(0);
353: }

355: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
356: {
358:   *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
359:   return(0);
360: }

362: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
363: {
365:   lim->ops->view    = PetscLimiterView_Sin;
366:   lim->ops->destroy = PetscLimiterDestroy_Sin;
367:   lim->ops->limit   = PetscLimiterLimit_Sin;
368:   return(0);
369: }

371: /*MC
372:   PETSCLIMITERSIN = "sin" - A PetscLimiter object

374:   Level: intermediate

376: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
377: M*/

379: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
380: {
381:   PetscLimiter_Sin *l;
382:   PetscErrorCode    ierr;

386:   PetscNewLog(lim, &l);
387:   lim->data = l;

389:   PetscLimiterInitialize_Sin(lim);
390:   return(0);
391: }

393: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
394: {
395:   PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
396:   PetscErrorCode    ierr;

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

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

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

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

422:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
423:   if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
424:   return(0);
425: }

427: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
428: {
430:   *phi = 0.0;
431:   return(0);
432: }

434: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
435: {
437:   lim->ops->view    = PetscLimiterView_Zero;
438:   lim->ops->destroy = PetscLimiterDestroy_Zero;
439:   lim->ops->limit   = PetscLimiterLimit_Zero;
440:   return(0);
441: }

443: /*MC
444:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

446:   Level: intermediate

448: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
449: M*/

451: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
452: {
453:   PetscLimiter_Zero *l;
454:   PetscErrorCode     ierr;

458:   PetscNewLog(lim, &l);
459:   lim->data = l;

461:   PetscLimiterInitialize_Zero(lim);
462:   return(0);
463: }

465: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
466: {
467:   PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
468:   PetscErrorCode    ierr;

471:   PetscFree(l);
472:   return(0);
473: }

475: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
476: {
477:   PetscViewerFormat format;
478:   PetscErrorCode    ierr;

481:   PetscViewerGetFormat(viewer, &format);
482:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
483:   return(0);
484: }

486: PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
487: {
488:   PetscBool      iascii;

494:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
495:   if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
496:   return(0);
497: }

499: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
500: {
502:   *phi = 1.0;
503:   return(0);
504: }

506: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
507: {
509:   lim->ops->view    = PetscLimiterView_None;
510:   lim->ops->destroy = PetscLimiterDestroy_None;
511:   lim->ops->limit   = PetscLimiterLimit_None;
512:   return(0);
513: }

515: /*MC
516:   PETSCLIMITERNONE = "none" - A PetscLimiter object

518:   Level: intermediate

520: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
521: M*/

523: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
524: {
525:   PetscLimiter_None *l;
526:   PetscErrorCode    ierr;

530:   PetscNewLog(lim, &l);
531:   lim->data = l;

533:   PetscLimiterInitialize_None(lim);
534:   return(0);
535: }

537: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
538: {
539:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
540:   PetscErrorCode    ierr;

543:   PetscFree(l);
544:   return(0);
545: }

547: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
548: {
549:   PetscViewerFormat format;
550:   PetscErrorCode    ierr;

553:   PetscViewerGetFormat(viewer, &format);
554:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
555:   return(0);
556: }

558: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
559: {
560:   PetscBool      iascii;

566:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
567:   if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
568:   return(0);
569: }

571: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
572: {
574:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
575:   return(0);
576: }

578: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
579: {
581:   lim->ops->view    = PetscLimiterView_Minmod;
582:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
583:   lim->ops->limit   = PetscLimiterLimit_Minmod;
584:   return(0);
585: }

587: /*MC
588:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

590:   Level: intermediate

592: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
593: M*/

595: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
596: {
597:   PetscLimiter_Minmod *l;
598:   PetscErrorCode    ierr;

602:   PetscNewLog(lim, &l);
603:   lim->data = l;

605:   PetscLimiterInitialize_Minmod(lim);
606:   return(0);
607: }

609: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
610: {
611:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
612:   PetscErrorCode    ierr;

615:   PetscFree(l);
616:   return(0);
617: }

619: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
620: {
621:   PetscViewerFormat format;
622:   PetscErrorCode    ierr;

625:   PetscViewerGetFormat(viewer, &format);
626:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
627:   return(0);
628: }

630: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
631: {
632:   PetscBool      iascii;

638:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
639:   if (iascii) {PetscLimiterView_VanLeer_Ascii(lim, viewer);}
640:   return(0);
641: }

643: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
644: {
646:   *phi = PetscMax(0, 4*f*(1-f));
647:   return(0);
648: }

650: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
651: {
653:   lim->ops->view    = PetscLimiterView_VanLeer;
654:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
655:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
656:   return(0);
657: }

659: /*MC
660:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

662:   Level: intermediate

664: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
665: M*/

667: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
668: {
669:   PetscLimiter_VanLeer *l;
670:   PetscErrorCode    ierr;

674:   PetscNewLog(lim, &l);
675:   lim->data = l;

677:   PetscLimiterInitialize_VanLeer(lim);
678:   return(0);
679: }

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

687:   PetscFree(l);
688:   return(0);
689: }

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

697:   PetscViewerGetFormat(viewer, &format);
698:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
699:   return(0);
700: }

702: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
703: {
704:   PetscBool      iascii;

710:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
711:   if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
712:   return(0);
713: }

715: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
716: {
718:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
719:   return(0);
720: }

722: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
723: {
725:   lim->ops->view    = PetscLimiterView_VanAlbada;
726:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
727:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
728:   return(0);
729: }

731: /*MC
732:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

734:   Level: intermediate

736: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
737: M*/

739: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
740: {
741:   PetscLimiter_VanAlbada *l;
742:   PetscErrorCode    ierr;

746:   PetscNewLog(lim, &l);
747:   lim->data = l;

749:   PetscLimiterInitialize_VanAlbada(lim);
750:   return(0);
751: }

753: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
754: {
755:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
756:   PetscErrorCode    ierr;

759:   PetscFree(l);
760:   return(0);
761: }

763: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
764: {
765:   PetscViewerFormat format;
766:   PetscErrorCode    ierr;

769:   PetscViewerGetFormat(viewer, &format);
770:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
771:   return(0);
772: }

774: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
775: {
776:   PetscBool      iascii;

782:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
783:   if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
784:   return(0);
785: }

787: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
788: {
790:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
791:   return(0);
792: }

794: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
795: {
797:   lim->ops->view    = PetscLimiterView_Superbee;
798:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
799:   lim->ops->limit   = PetscLimiterLimit_Superbee;
800:   return(0);
801: }

803: /*MC
804:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

806:   Level: intermediate

808: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
809: M*/

811: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
812: {
813:   PetscLimiter_Superbee *l;
814:   PetscErrorCode    ierr;

818:   PetscNewLog(lim, &l);
819:   lim->data = l;

821:   PetscLimiterInitialize_Superbee(lim);
822:   return(0);
823: }

825: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
826: {
827:   PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
828:   PetscErrorCode    ierr;

831:   PetscFree(l);
832:   return(0);
833: }

835: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
836: {
837:   PetscViewerFormat format;
838:   PetscErrorCode    ierr;

841:   PetscViewerGetFormat(viewer, &format);
842:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
843:   return(0);
844: }

846: PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
847: {
848:   PetscBool      iascii;

854:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
855:   if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
856:   return(0);
857: }

859: /* aka Barth-Jespersen */
860: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
861: {
863:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
864:   return(0);
865: }

867: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
868: {
870:   lim->ops->view    = PetscLimiterView_MC;
871:   lim->ops->destroy = PetscLimiterDestroy_MC;
872:   lim->ops->limit   = PetscLimiterLimit_MC;
873:   return(0);
874: }

876: /*MC
877:   PETSCLIMITERMC = "mc" - A PetscLimiter object

879:   Level: intermediate

881: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
882: M*/

884: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
885: {
886:   PetscLimiter_MC *l;
887:   PetscErrorCode    ierr;

891:   PetscNewLog(lim, &l);
892:   lim->data = l;

894:   PetscLimiterInitialize_MC(lim);
895:   return(0);
896: }

898: PetscClassId PETSCFV_CLASSID = 0;

900: PetscFunctionList PetscFVList              = NULL;
901: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

903: /*@C
904:   PetscFVRegister - Adds a new PetscFV implementation

906:   Not Collective

908:   Input Parameters:
909: + name        - The name of a new user-defined creation routine
910: - create_func - The creation routine itself

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

915:   Sample usage:
916: .vb
917:     PetscFVRegister("my_fv", MyPetscFVCreate);
918: .ve

920:   Then, your PetscFV type can be chosen with the procedural interface via
921: .vb
922:     PetscFVCreate(MPI_Comm, PetscFV *);
923:     PetscFVSetType(PetscFV, "my_fv");
924: .ve
925:    or at runtime via the option
926: .vb
927:     -petscfv_type my_fv
928: .ve

930:   Level: advanced

932: .keywords: PetscFV, register
933: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()

935: @*/
936: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
937: {

941:   PetscFunctionListAdd(&PetscFVList, sname, function);
942:   return(0);
943: }

945: /*@C
946:   PetscFVSetType - Builds a particular PetscFV

948:   Collective on PetscFV

950:   Input Parameters:
951: + fvm  - The PetscFV object
952: - name - The kind of FVM space

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

957:   Level: intermediate

959: .keywords: PetscFV, set, type
960: .seealso: PetscFVGetType(), PetscFVCreate()
961: @*/
962: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
963: {
964:   PetscErrorCode (*r)(PetscFV);
965:   PetscBool      match;

970:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
971:   if (match) return(0);

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

977:   if (fvm->ops->destroy) {
978:     (*fvm->ops->destroy)(fvm);
979:     fvm->ops->destroy = NULL;
980:   }
981:   (*r)(fvm);
982:   PetscObjectChangeTypeName((PetscObject) fvm, name);
983:   return(0);
984: }

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

989:   Not Collective

991:   Input Parameter:
992: . fvm  - The PetscFV

994:   Output Parameter:
995: . name - The PetscFV type name

997:   Level: intermediate

999: .keywords: PetscFV, get, type, name
1000: .seealso: PetscFVSetType(), PetscFVCreate()
1001: @*/
1002: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1003: {

1009:   PetscFVRegisterAll();
1010:   *name = ((PetscObject) fvm)->type_name;
1011:   return(0);
1012: }

1014: /*@C
1015:   PetscFVView - Views a PetscFV

1017:   Collective on PetscFV

1019:   Input Parameter:
1020: + fvm - the PetscFV object to view
1021: - v   - the viewer

1023:   Level: developer

1025: .seealso: PetscFVDestroy()
1026: @*/
1027: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1028: {

1033:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1034:   if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1035:   return(0);
1036: }

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

1041:   Collective on PetscFV

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

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

1049:   Level: developer

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

1062:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1063:   else                                 defaultType = ((PetscObject) fvm)->type_name;
1064:   PetscFVRegisterAll();

1066:   PetscObjectOptionsBegin((PetscObject) fvm);
1067:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1068:   if (flg) {
1069:     PetscFVSetType(fvm, name);
1070:   } else if (!((PetscObject) fvm)->type_name) {
1071:     PetscFVSetType(fvm, defaultType);

1073:   }
1074:   PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1075:   if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1076:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1077:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1078:   PetscLimiterSetFromOptions(fvm->limiter);
1079:   PetscOptionsEnd();
1080:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1081:   return(0);
1082: }

1084: /*@
1085:   PetscFVSetUp - Construct data structures for the PetscFV

1087:   Collective on PetscFV

1089:   Input Parameter:
1090: . fvm - the PetscFV object to setup

1092:   Level: developer

1094: .seealso: PetscFVView(), PetscFVDestroy()
1095: @*/
1096: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1097: {

1102:   PetscLimiterSetUp(fvm->limiter);
1103:   if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1104:   return(0);
1105: }

1107: /*@
1108:   PetscFVDestroy - Destroys a PetscFV object

1110:   Collective on PetscFV

1112:   Input Parameter:
1113: . fvm - the PetscFV object to destroy

1115:   Level: developer

1117: .seealso: PetscFVView()
1118: @*/
1119: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1120: {
1121:   PetscInt       i;

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

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

1131:   for (i = 0; i < (*fvm)->numComponents; i++) {
1132:     PetscFree((*fvm)->componentNames[i]);
1133:   }
1134:   PetscFree((*fvm)->componentNames);
1135:   PetscLimiterDestroy(&(*fvm)->limiter);
1136:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1137:   PetscFree((*fvm)->fluxWork);
1138:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1139:   PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);

1141:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1142:   PetscHeaderDestroy(fvm);
1143:   return(0);
1144: }

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

1149:   Collective on MPI_Comm

1151:   Input Parameter:
1152: . comm - The communicator for the PetscFV object

1154:   Output Parameter:
1155: . fvm - The PetscFV object

1157:   Level: beginner

1159: .seealso: PetscFVSetType(), PETSCFVUPWIND
1160: @*/
1161: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1162: {
1163:   PetscFV        f;

1168:   *fvm = NULL;
1169:   PetscFVInitializePackage();

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

1174:   PetscLimiterCreate(comm, &f->limiter);
1175:   f->numComponents    = 1;
1176:   f->dim              = 0;
1177:   f->computeGradients = PETSC_FALSE;
1178:   f->fluxWork         = NULL;
1179:   PetscCalloc1(f->numComponents,&f->componentNames);

1181:   *fvm = f;
1182:   return(0);
1183: }

1185: /*@
1186:   PetscFVSetLimiter - Set the limiter object

1188:   Logically collective on PetscFV

1190:   Input Parameters:
1191: + fvm - the PetscFV object
1192: - lim - The PetscLimiter

1194:   Level: developer

1196: .seealso: PetscFVGetLimiter()
1197: @*/
1198: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1199: {

1205:   PetscLimiterDestroy(&fvm->limiter);
1206:   PetscObjectReference((PetscObject) lim);
1207:   fvm->limiter = lim;
1208:   return(0);
1209: }

1211: /*@
1212:   PetscFVGetLimiter - Get the limiter object

1214:   Not collective

1216:   Input Parameter:
1217: . fvm - the PetscFV object

1219:   Output Parameter:
1220: . lim - The PetscLimiter

1222:   Level: developer

1224: .seealso: PetscFVSetLimiter()
1225: @*/
1226: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1227: {
1231:   *lim = fvm->limiter;
1232:   return(0);
1233: }

1235: /*@
1236:   PetscFVSetNumComponents - Set the number of field components

1238:   Logically collective on PetscFV

1240:   Input Parameters:
1241: + fvm - the PetscFV object
1242: - comp - The number of components

1244:   Level: developer

1246: .seealso: PetscFVGetNumComponents()
1247: @*/
1248: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1249: {

1254:   if (fvm->numComponents != comp) {
1255:     PetscInt i;

1257:     for (i = 0; i < fvm->numComponents; i++) {
1258:       PetscFree(fvm->componentNames[i]);
1259:     }
1260:     PetscFree(fvm->componentNames);
1261:     PetscCalloc1(comp,&fvm->componentNames);
1262:   }
1263:   fvm->numComponents = comp;
1264:   PetscFree(fvm->fluxWork);
1265:   PetscMalloc1(comp, &fvm->fluxWork);
1266:   return(0);
1267: }

1269: /*@
1270:   PetscFVGetNumComponents - Get the number of field components

1272:   Not collective

1274:   Input Parameter:
1275: . fvm - the PetscFV object

1277:   Output Parameter:
1278: , comp - The number of components

1280:   Level: developer

1282: .seealso: PetscFVSetNumComponents()
1283: @*/
1284: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1285: {
1289:   *comp = fvm->numComponents;
1290:   return(0);
1291: }

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

1296:   Logically collective on PetscFV
1297:   Input Parameters:
1298: + fvm - the PetscFV object
1299: . comp - the component number
1300: - name - the component name

1302:   Level: developer

1304: .seealso: PetscFVGetComponentName()
1305: @*/
1306: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1307: {

1311:   PetscFree(fvm->componentNames[comp]);
1312:   PetscStrallocpy(name,&fvm->componentNames[comp]);
1313:   return(0);
1314: }

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

1319:   Logically collective on PetscFV
1320:   Input Parameters:
1321: + fvm - the PetscFV object
1322: - comp - the component number

1324:   Output Parameter:
1325: . name - the component name

1327:   Level: developer

1329: .seealso: PetscFVSetComponentName()
1330: @*/
1331: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1332: {
1334:   *name = fvm->componentNames[comp];
1335:   return(0);
1336: }

1338: /*@
1339:   PetscFVSetSpatialDimension - Set the spatial dimension

1341:   Logically collective on PetscFV

1343:   Input Parameters:
1344: + fvm - the PetscFV object
1345: - dim - The spatial dimension

1347:   Level: developer

1349: .seealso: PetscFVGetSpatialDimension()
1350: @*/
1351: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1352: {
1355:   fvm->dim = dim;
1356:   return(0);
1357: }

1359: /*@
1360:   PetscFVGetSpatialDimension - Get the spatial dimension

1362:   Logically collective on PetscFV

1364:   Input Parameter:
1365: . fvm - the PetscFV object

1367:   Output Parameter:
1368: . dim - The spatial dimension

1370:   Level: developer

1372: .seealso: PetscFVSetSpatialDimension()
1373: @*/
1374: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1375: {
1379:   *dim = fvm->dim;
1380:   return(0);
1381: }

1383: /*@
1384:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1386:   Logically collective on PetscFV

1388:   Input Parameters:
1389: + fvm - the PetscFV object
1390: - computeGradients - Flag to compute cell gradients

1392:   Level: developer

1394: .seealso: PetscFVGetComputeGradients()
1395: @*/
1396: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1397: {
1400:   fvm->computeGradients = computeGradients;
1401:   return(0);
1402: }

1404: /*@
1405:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1407:   Not collective

1409:   Input Parameter:
1410: . fvm - the PetscFV object

1412:   Output Parameter:
1413: . computeGradients - Flag to compute cell gradients

1415:   Level: developer

1417: .seealso: PetscFVSetComputeGradients()
1418: @*/
1419: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1420: {
1424:   *computeGradients = fvm->computeGradients;
1425:   return(0);
1426: }

1428: /*@
1429:   PetscFVSetQuadrature - Set the quadrature object

1431:   Logically collective on PetscFV

1433:   Input Parameters:
1434: + fvm - the PetscFV object
1435: - q - The PetscQuadrature

1437:   Level: developer

1439: .seealso: PetscFVGetQuadrature()
1440: @*/
1441: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1442: {

1447:   PetscQuadratureDestroy(&fvm->quadrature);
1448:   PetscObjectReference((PetscObject) q);
1449:   fvm->quadrature = q;
1450:   return(0);
1451: }

1453: /*@
1454:   PetscFVGetQuadrature - Get the quadrature object

1456:   Not collective

1458:   Input Parameter:
1459: . fvm - the PetscFV object

1461:   Output Parameter:
1462: . lim - The PetscQuadrature

1464:   Level: developer

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

1478:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1479:     PetscCalloc1(fvm->dim, &points);
1480:     PetscMalloc1(1, &weights);
1481:     weights[0] = 1.0;
1482:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1483:   }
1484:   *q = fvm->quadrature;
1485:   return(0);
1486: }

1488: /*@
1489:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1491:   Not collective

1493:   Input Parameter:
1494: . fvm - The PetscFV object

1496:   Output Parameter:
1497: . sp - The PetscDualSpace object

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

1501:   Level: developer

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

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

1530:       PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1531:       PetscCalloc1(dim, &points);
1532:       PetscCalloc1(Nc, &weights);
1533:       weights[c] = 1.0;
1534:       PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1535:       PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1536:       PetscQuadratureDestroy(&qc);
1537:     }
1538:   }
1539:   *sp = fvm->dualSpace;
1540:   return(0);
1541: }

1543: /*@
1544:   PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product

1546:   Not collective

1548:   Input Parameters:
1549: + fvm - The PetscFV object
1550: - sp  - The PetscDualSpace object

1552:   Level: intermediate

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

1556: .seealso: PetscFVCreate()
1557: @*/
1558: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1559: {

1565:   PetscDualSpaceDestroy(&fvm->dualSpace);
1566:   fvm->dualSpace = sp;
1567:   PetscObjectReference((PetscObject) fvm->dualSpace);
1568:   return(0);
1569: }

1571: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1572: {
1573:   PetscInt         npoints;
1574:   const PetscReal *points;
1575:   PetscErrorCode   ierr;

1582:   PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1583:   if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1584:   if (B) *B = fvm->B;
1585:   if (D) *D = fvm->D;
1586:   if (H) *H = fvm->H;
1587:   return(0);
1588: }

1590: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1591: {
1592:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1593:   PetscInt         dim;      /* Spatial dimension */
1594:   PetscInt         comp;     /* Field components */
1595:   PetscInt         p, d, c, e;
1596:   PetscErrorCode   ierr;

1604:   PetscFVGetSpatialDimension(fvm, &dim);
1605:   PetscFVGetNumComponents(fvm, &comp);
1606:   if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1607:   if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1608:   if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1609:   if (B) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) (*B)[(p*pdim + d)*comp + c] = 1.0;}
1610:   if (D) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim; ++e) (*D)[((p*pdim + d)*comp + c)*dim + e] = 0.0;}
1611:   if (H) {for (p = 0; p < npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < comp; ++c) for (e = 0; e < dim*dim; ++e) (*H)[((p*pdim + d)*comp + c)*dim*dim + e] = 0.0;}
1612:   return(0);
1613: }

1615: PetscErrorCode PetscFVRestoreTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1616: {

1621:   if (B && *B) {PetscFree(*B);}
1622:   if (D && *D) {PetscFree(*D);}
1623:   if (H && *H) {PetscFree(*H);}
1624:   return(0);
1625: }

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

1630:   Input Parameters:
1631: + fvm      - The PetscFV object
1632: . numFaces - The number of cell faces which are not constrained
1633: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1635:   Level: developer

1637: .seealso: PetscFVCreate()
1638: @*/
1639: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1640: {

1645:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1646:   return(0);
1647: }

1649: /*C
1650:   PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration

1652:   Not collective

1654:   Input Parameters:
1655: + fvm          - The PetscFV object for the field being integrated
1656: . prob         - The PetscDS specifing the discretizations and continuum functions
1657: . field        - The field being integrated
1658: . Nf           - The number of faces in the chunk
1659: . fgeom        - The face geometry for each face in the chunk
1660: . neighborVol  - The volume for each pair of cells in the chunk
1661: . uL           - The state from the cell on the left
1662: - uR           - The state from the cell on the right

1664:   Output Parameter
1665: + fluxL        - the left fluxes for each face
1666: - fluxR        - the right fluxes for each face
1667: */
1668: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1669:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1670: {

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

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

1684:   Input Parameter:
1685: . fv - The initial PetscFV

1687:   Output Parameter:
1688: . fvRef - The refined PetscFV

1690:   Level: developer

1692: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1693: @*/
1694: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1695: {
1696:   PetscDualSpace   Q, Qref;
1697:   DM               K, Kref;
1698:   PetscQuadrature  q, qref;
1699:   CellRefiner      cellRefiner;
1700:   PetscReal       *v0;
1701:   PetscReal       *jac, *invjac;
1702:   PetscInt         numComp, numSubelements, s;
1703:   PetscErrorCode   ierr;

1706:   PetscFVGetDualSpace(fv, &Q);
1707:   PetscFVGetQuadrature(fv, &q);
1708:   PetscDualSpaceGetDM(Q, &K);
1709:   /* Create dual space */
1710:   PetscDualSpaceDuplicate(Q, &Qref);
1711:   DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1712:   PetscDualSpaceSetDM(Qref, Kref);
1713:   DMDestroy(&Kref);
1714:   PetscDualSpaceSetUp(Qref);
1715:   /* Create volume */
1716:   PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1717:   PetscFVSetDualSpace(*fvRef, Qref);
1718:   PetscFVGetNumComponents(fv,    &numComp);
1719:   PetscFVSetNumComponents(*fvRef, numComp);
1720:   PetscFVSetUp(*fvRef);
1721:   /* Create quadrature */
1722:   DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1723:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1724:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1725:   PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1726:   for (s = 0; s < numSubelements; ++s) {
1727:     PetscQuadrature  qs;
1728:     const PetscReal *points, *weights;
1729:     PetscReal       *p, *w;
1730:     PetscInt         dim, Nc, npoints, np;

1732:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1733:     PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1734:     np   = npoints/numSubelements;
1735:     PetscMalloc1(np*dim,&p);
1736:     PetscMalloc1(np*Nc,&w);
1737:     PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1738:     PetscMemcpy(w, &weights[s*np*Nc], np*Nc  * sizeof(PetscReal));
1739:     PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1740:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1741:     PetscQuadratureDestroy(&qs);
1742:   }
1743:   CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1744:   PetscFVSetQuadrature(*fvRef, qref);
1745:   PetscQuadratureDestroy(&qref);
1746:   PetscDualSpaceDestroy(&Qref);
1747:   return(0);
1748: }

1750: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1751: {
1752:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1756:   PetscFree(b);
1757:   return(0);
1758: }

1760: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1761: {
1762:   PetscInt          Nc, c;
1763:   PetscViewerFormat format;
1764:   PetscErrorCode    ierr;

1767:   PetscFVGetNumComponents(fv, &Nc);
1768:   PetscViewerGetFormat(viewer, &format);
1769:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1770:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1771:   for (c = 0; c < Nc; c++) {
1772:     if (fv->componentNames[c]) {
1773:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1774:     }
1775:   }
1776:   return(0);
1777: }

1779: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1780: {
1781:   PetscBool      iascii;

1787:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1788:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1789:   return(0);
1790: }

1792: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1793: {
1795:   return(0);
1796: }

1798: /*
1799:   neighborVol[f*2+0] contains the left  geom
1800:   neighborVol[f*2+1] contains the right geom
1801: */
1802: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1803:                                                   PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1804: {
1805:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1806:   void              *rctx;
1807:   PetscScalar       *flux = fvm->fluxWork;
1808:   const PetscScalar *constants;
1809:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;
1810:   PetscErrorCode     ierr;

1813:   PetscDSGetTotalComponents(prob, &Nc);
1814:   PetscDSGetTotalDimension(prob, &totDim);
1815:   PetscDSGetFieldOffset(prob, field, &off);
1816:   PetscDSGetRiemannSolver(prob, field, &riemann);
1817:   PetscDSGetContext(prob, field, &rctx);
1818:   PetscDSGetConstants(prob, &numConstants, &constants);
1819:   PetscFVGetSpatialDimension(fvm, &dim);
1820:   PetscFVGetNumComponents(fvm, &pdim);
1821:   for (f = 0; f < Nf; ++f) {
1822:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1823:     for (d = 0; d < pdim; ++d) {
1824:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1825:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1826:     }
1827:   }
1828:   return(0);
1829: }

1831: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1832: {
1834:   fvm->ops->setfromoptions          = NULL;
1835:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1836:   fvm->ops->view                    = PetscFVView_Upwind;
1837:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
1838:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
1839:   return(0);
1840: }

1842: /*MC
1843:   PETSCFVUPWIND = "upwind" - A PetscFV object

1845:   Level: intermediate

1847: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1848: M*/

1850: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1851: {
1852:   PetscFV_Upwind *b;
1853:   PetscErrorCode  ierr;

1857:   PetscNewLog(fvm,&b);
1858:   fvm->data = b;

1860:   PetscFVInitialize_Upwind(fvm);
1861:   return(0);
1862: }

1864:  #include <petscblaslapack.h>

1866: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1867: {
1868:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1869:   PetscErrorCode        ierr;

1872:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1873:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1874:   PetscFree(ls);
1875:   return(0);
1876: }

1878: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1879: {
1880:   PetscInt          Nc, c;
1881:   PetscViewerFormat format;
1882:   PetscErrorCode    ierr;

1885:   PetscFVGetNumComponents(fv, &Nc);
1886:   PetscViewerGetFormat(viewer, &format);
1887:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1888:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1889:   for (c = 0; c < Nc; c++) {
1890:     if (fv->componentNames[c]) {
1891:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1892:     }
1893:   }
1894:   return(0);
1895: }

1897: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1898: {
1899:   PetscBool      iascii;

1905:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1906:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1907:   return(0);
1908: }

1910: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1911: {
1913:   return(0);
1914: }

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

1925:   if (debug) {
1926:     PetscMalloc1(m*n,&Aback);
1927:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1928:   }

1930:   PetscBLASIntCast(m,&M);
1931:   PetscBLASIntCast(n,&N);
1932:   PetscBLASIntCast(mstride,&lda);
1933:   PetscBLASIntCast(worksize,&ldwork);
1934:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1935:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1936:   PetscFPTrapPop();
1937:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1938:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1940:   /* Extract an explicit representation of Q */
1941:   Q    = Ainv;
1942:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1943:   K    = N;                     /* full rank */
1944:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1945:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1947:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1948:   Alpha = 1.0;
1949:   ldb   = lda;
1950:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1951:   /* Ainv is Q, overwritten with inverse */

1953:   if (debug) {                      /* Check that pseudo-inverse worked */
1954:     PetscScalar  Beta = 0.0;
1955:     PetscBLASInt ldc;
1956:     K   = N;
1957:     ldc = N;
1958:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1959:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1960:     PetscFree(Aback);
1961:   }
1962:   return(0);
1963: }

1965: /* Overwrites A. Can handle degenerate problems and m<n. */
1966: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1967: {
1968:   PetscBool      debug = PETSC_FALSE;
1969:   PetscScalar   *Brhs, *Aback;
1970:   PetscScalar   *tmpwork;
1971:   PetscReal      rcond;
1972: #if defined (PETSC_USE_COMPLEX)
1973:   PetscInt       rworkSize;
1974:   PetscReal     *rwork;
1975: #endif
1976:   PetscInt       i, j, maxmn;
1977:   PetscBLASInt   M, N, lda, ldb, ldwork;
1978:   PetscBLASInt   nrhs, irank, info;

1982:   if (debug) {
1983:     PetscMalloc1(m*n,&Aback);
1984:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1985:   }

1987:   /* initialize to identity */
1988:   tmpwork = Ainv;
1989:   Brhs = work;
1990:   maxmn = PetscMax(m,n);
1991:   for (j=0; j<maxmn; j++) {
1992:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1993:   }

1995:   PetscBLASIntCast(m,&M);
1996:   PetscBLASIntCast(n,&N);
1997:   PetscBLASIntCast(mstride,&lda);
1998:   PetscBLASIntCast(maxmn,&ldb);
1999:   PetscBLASIntCast(worksize,&ldwork);
2000:   rcond = -1;
2001:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2002:   nrhs  = M;
2003: #if defined(PETSC_USE_COMPLEX)
2004:   rworkSize = 5 * PetscMin(M,N);
2005:   PetscMalloc1(rworkSize,&rwork);
2006:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2007:   PetscFPTrapPop();
2008:   PetscFree(rwork);
2009: #else
2010:   nrhs  = M;
2011:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2012:   PetscFPTrapPop();
2013: #endif
2014:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2015:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2016:   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");

2018:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2019:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2020:   for (i=0; i<n; i++) {
2021:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2022:   }
2023:   return(0);
2024: }

2026: #if 0
2027: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2028: {
2029:   PetscReal       grad[2] = {0, 0};
2030:   const PetscInt *faces;
2031:   PetscInt        numFaces, f;
2032:   PetscErrorCode  ierr;

2035:   DMPlexGetConeSize(dm, cell, &numFaces);
2036:   DMPlexGetCone(dm, cell, &faces);
2037:   for (f = 0; f < numFaces; ++f) {
2038:     const PetscInt *fcells;
2039:     const CellGeom *cg1;
2040:     const FaceGeom *fg;

2042:     DMPlexGetSupport(dm, faces[f], &fcells);
2043:     DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2044:     for (i = 0; i < 2; ++i) {
2045:       PetscScalar du;

2047:       if (fcells[i] == c) continue;
2048:       DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2049:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2050:       grad[0] += fg->grad[!i][0] * du;
2051:       grad[1] += fg->grad[!i][1] * du;
2052:     }
2053:   }
2054:   PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2055:   return(0);
2056: }
2057: #endif

2059: /*
2060:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

2062:   Input Parameters:
2063: + fvm      - The PetscFV object
2064: . numFaces - The number of cell faces which are not constrained
2065: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

2067:   Level: developer

2069: .seealso: PetscFVCreate()
2070: */
2071: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2072: {
2073:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *) fvm->data;
2074:   const PetscBool       useSVD   = PETSC_TRUE;
2075:   const PetscInt        maxFaces = ls->maxFaces;
2076:   PetscInt              dim, f, d;
2077:   PetscErrorCode        ierr;

2080:   if (numFaces > maxFaces) {
2081:     if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2082:     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2083:   }
2084:   PetscFVGetSpatialDimension(fvm, &dim);
2085:   for (f = 0; f < numFaces; ++f) {
2086:     for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2087:   }
2088:   /* Overwrites B with garbage, returns Binv in row-major format */
2089:   if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2090:   else        {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2091:   for (f = 0; f < numFaces; ++f) {
2092:     for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2093:   }
2094:   return(0);
2095: }

2097: /*
2098:   neighborVol[f*2+0] contains the left  geom
2099:   neighborVol[f*2+1] contains the right geom
2100: */
2101: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2102:                                                         PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2103: {
2104:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2105:   void              *rctx;
2106:   PetscScalar       *flux = fvm->fluxWork;
2107:   const PetscScalar *constants;
2108:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;
2109:   PetscErrorCode     ierr;

2112:   PetscDSGetTotalComponents(prob, &Nc);
2113:   PetscDSGetTotalDimension(prob, &totDim);
2114:   PetscDSGetFieldOffset(prob, field, &off);
2115:   PetscDSGetRiemannSolver(prob, field, &riemann);
2116:   PetscDSGetContext(prob, field, &rctx);
2117:   PetscDSGetConstants(prob, &numConstants, &constants);
2118:   PetscFVGetSpatialDimension(fvm, &dim);
2119:   PetscFVGetNumComponents(fvm, &pdim);
2120:   for (f = 0; f < Nf; ++f) {
2121:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2122:     for (d = 0; d < pdim; ++d) {
2123:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2124:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2125:     }
2126:   }
2127:   return(0);
2128: }

2130: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2131: {
2132:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2133:   PetscInt              dim, m, n, nrhs, minwork;
2134:   PetscErrorCode        ierr;

2138:   PetscFVGetSpatialDimension(fvm, &dim);
2139:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2140:   ls->maxFaces = maxFaces;
2141:   m       = ls->maxFaces;
2142:   n       = dim;
2143:   nrhs    = ls->maxFaces;
2144:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2145:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2146:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2147:   return(0);
2148: }

2150: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2151: {
2153:   fvm->ops->setfromoptions          = NULL;
2154:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2155:   fvm->ops->view                    = PetscFVView_LeastSquares;
2156:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2157:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2158:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2159:   return(0);
2160: }

2162: /*MC
2163:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2165:   Level: intermediate

2167: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2168: M*/

2170: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2171: {
2172:   PetscFV_LeastSquares *ls;
2173:   PetscErrorCode        ierr;

2177:   PetscNewLog(fvm, &ls);
2178:   fvm->data = ls;

2180:   ls->maxFaces = -1;
2181:   ls->workSize = -1;
2182:   ls->B        = NULL;
2183:   ls->Binv     = NULL;
2184:   ls->tau      = NULL;
2185:   ls->work     = NULL;

2187:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2188:   PetscFVInitialize_LeastSquares(fvm);
2189:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2190:   return(0);
2191: }

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

2196:   Not collective

2198:   Input parameters:
2199: + fvm      - The PetscFV object
2200: - maxFaces - The maximum number of cell faces

2202:   Level: intermediate

2204: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2205: @*/
2206: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2207: {

2212:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2213:   return(0);
2214: }