Actual source code: dtfv.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/petscfvimpl.h> /*I "petscfv.h" I*/
  2: #include <petsc/private/dmpleximpl.h> /* For CellRefiner */
  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";

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

 23:   Not Collective

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

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

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

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

 47:   Level: advanced

 49: .keywords: PetscLimiter, register
 50: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()

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

 58:   PetscFunctionListAdd(&PetscLimiterList, sname, function);
 59:   return(0);
 60: }

 64: /*@C
 65:   PetscLimiterSetType - Builds a particular PetscLimiter

 67:   Collective on PetscLimiter

 69:   Input Parameters:
 70: + lim  - The PetscLimiter object
 71: - name - The kind of limiter

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

 76:   Level: intermediate

 78: .keywords: PetscLimiter, set, type
 79: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
 80: @*/
 81: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
 82: {
 83:   PetscErrorCode (*r)(PetscLimiter);
 84:   PetscBool      match;

 89:   PetscObjectTypeCompare((PetscObject) lim, name, &match);
 90:   if (match) return(0);

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

 96:   if (lim->ops->destroy) {
 97:     (*lim->ops->destroy)(lim);
 98:     lim->ops->destroy = NULL;
 99:   }
100:   (*r)(lim);
101:   PetscObjectChangeTypeName((PetscObject) lim, name);
102:   return(0);
103: }

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

110:   Not Collective

112:   Input Parameter:
113: . lim  - The PetscLimiter

115:   Output Parameter:
116: . name - The PetscLimiter type name

118:   Level: intermediate

120: .keywords: PetscLimiter, get, type, name
121: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
122: @*/
123: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
124: {

130:   PetscLimiterRegisterAll();
131:   *name = ((PetscObject) lim)->type_name;
132:   return(0);
133: }

137: /*@C
138:   PetscLimiterView - Views a PetscLimiter

140:   Collective on PetscLimiter

142:   Input Parameter:
143: + lim - the PetscLimiter object to view
144: - v   - the viewer

146:   Level: developer

148: .seealso: PetscLimiterDestroy()
149: @*/
150: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
151: {

156:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);}
157:   if (lim->ops->view) {(*lim->ops->view)(lim, v);}
158:   return(0);
159: }

163: /*@
164:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

166:   Collective on PetscLimiter

168:   Input Parameter:
169: . lim - the PetscLimiter object to set options for

171:   Level: developer

173: .seealso: PetscLimiterView()
174: @*/
175: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
176: {
177:   const char    *defaultType;
178:   char           name[256];
179:   PetscBool      flg;

184:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
185:   else                                 defaultType = ((PetscObject) lim)->type_name;
186:   PetscLimiterRegisterAll();

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

205: /*@C
206:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

208:   Collective on PetscLimiter

210:   Input Parameter:
211: . lim - the PetscLimiter object to setup

213:   Level: developer

215: .seealso: PetscLimiterView(), PetscLimiterDestroy()
216: @*/
217: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
218: {

223:   if (lim->ops->setup) {(*lim->ops->setup)(lim);}
224:   return(0);
225: }

229: /*@
230:   PetscLimiterDestroy - Destroys a PetscLimiter object

232:   Collective on PetscLimiter

234:   Input Parameter:
235: . lim - the PetscLimiter object to destroy

237:   Level: developer

239: .seealso: PetscLimiterView()
240: @*/
241: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
242: {

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

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

252:   if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
253:   PetscHeaderDestroy(lim);
254:   return(0);
255: }

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

262:   Collective on MPI_Comm

264:   Input Parameter:
265: . comm - The communicator for the PetscLimiter object

267:   Output Parameter:
268: . lim - The PetscLimiter object

270:   Level: beginner

272: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
273: @*/
274: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
275: {
276:   PetscLimiter   l;

281:   PetscCitationsRegister(LimiterCitation,&Limitercite);
282:   *lim = NULL;
283:   PetscFVInitializePackage();

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

287:   *lim = l;
288:   return(0);
289: }

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

335:   (*lim->ops->limit)(lim, flim, phi);
336:   return(0);
337: }

341: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
342: {
343:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
344:   PetscErrorCode    ierr;

347:   PetscFree(l);
348:   return(0);
349: }

353: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
354: {
355:   PetscViewerFormat format;
356:   PetscErrorCode    ierr;

359:   PetscViewerGetFormat(viewer, &format);
360:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
361:   return(0);
362: }

366: PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
367: {
368:   PetscBool      iascii;

374:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
375:   if (iascii) {PetscLimiterView_Sin_Ascii(lim, viewer);}
376:   return(0);
377: }

381: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
382: {
384:   *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
385:   return(0);
386: }

390: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
391: {
393:   lim->ops->view    = PetscLimiterView_Sin;
394:   lim->ops->destroy = PetscLimiterDestroy_Sin;
395:   lim->ops->limit   = PetscLimiterLimit_Sin;
396:   return(0);
397: }

399: /*MC
400:   PETSCLIMITERSIN = "sin" - A PetscLimiter object

402:   Level: intermediate

404: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
405: M*/

409: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
410: {
411:   PetscLimiter_Sin *l;
412:   PetscErrorCode    ierr;

416:   PetscNewLog(lim, &l);
417:   lim->data = l;

419:   PetscLimiterInitialize_Sin(lim);
420:   return(0);
421: }

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

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

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

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

450: PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
451: {
452:   PetscBool      iascii;

458:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
459:   if (iascii) {PetscLimiterView_Zero_Ascii(lim, viewer);}
460:   return(0);
461: }

465: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
466: {
468:   *phi = 0.0;
469:   return(0);
470: }

474: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
475: {
477:   lim->ops->view    = PetscLimiterView_Zero;
478:   lim->ops->destroy = PetscLimiterDestroy_Zero;
479:   lim->ops->limit   = PetscLimiterLimit_Zero;
480:   return(0);
481: }

483: /*MC
484:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

486:   Level: intermediate

488: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
489: M*/

493: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
494: {
495:   PetscLimiter_Zero *l;
496:   PetscErrorCode     ierr;

500:   PetscNewLog(lim, &l);
501:   lim->data = l;

503:   PetscLimiterInitialize_Zero(lim);
504:   return(0);
505: }

509: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
510: {
511:   PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
512:   PetscErrorCode    ierr;

515:   PetscFree(l);
516:   return(0);
517: }

521: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
522: {
523:   PetscViewerFormat format;
524:   PetscErrorCode    ierr;

527:   PetscViewerGetFormat(viewer, &format);
528:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
529:   return(0);
530: }

534: PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
535: {
536:   PetscBool      iascii;

542:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
543:   if (iascii) {PetscLimiterView_None_Ascii(lim, viewer);}
544:   return(0);
545: }

549: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
550: {
552:   *phi = 1.0;
553:   return(0);
554: }

558: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
559: {
561:   lim->ops->view    = PetscLimiterView_None;
562:   lim->ops->destroy = PetscLimiterDestroy_None;
563:   lim->ops->limit   = PetscLimiterLimit_None;
564:   return(0);
565: }

567: /*MC
568:   PETSCLIMITERNONE = "none" - A PetscLimiter object

570:   Level: intermediate

572: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
573: M*/

577: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
578: {
579:   PetscLimiter_None *l;
580:   PetscErrorCode    ierr;

584:   PetscNewLog(lim, &l);
585:   lim->data = l;

587:   PetscLimiterInitialize_None(lim);
588:   return(0);
589: }

593: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
594: {
595:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
596:   PetscErrorCode    ierr;

599:   PetscFree(l);
600:   return(0);
601: }

605: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
606: {
607:   PetscViewerFormat format;
608:   PetscErrorCode    ierr;

611:   PetscViewerGetFormat(viewer, &format);
612:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
613:   return(0);
614: }

618: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
619: {
620:   PetscBool      iascii;

626:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
627:   if (iascii) {PetscLimiterView_Minmod_Ascii(lim, viewer);}
628:   return(0);
629: }

633: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
634: {
636:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
637:   return(0);
638: }

642: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
643: {
645:   lim->ops->view    = PetscLimiterView_Minmod;
646:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
647:   lim->ops->limit   = PetscLimiterLimit_Minmod;
648:   return(0);
649: }

651: /*MC
652:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

654:   Level: intermediate

656: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
657: M*/

661: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
662: {
663:   PetscLimiter_Minmod *l;
664:   PetscErrorCode    ierr;

668:   PetscNewLog(lim, &l);
669:   lim->data = l;

671:   PetscLimiterInitialize_Minmod(lim);
672:   return(0);
673: }

677: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
678: {
679:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
680:   PetscErrorCode    ierr;

683:   PetscFree(l);
684:   return(0);
685: }

689: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
690: {
691:   PetscViewerFormat format;
692:   PetscErrorCode    ierr;

695:   PetscViewerGetFormat(viewer, &format);
696:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
697:   return(0);
698: }

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

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

717: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
718: {
720:   *phi = PetscMax(0, 4*f*(1-f));
721:   return(0);
722: }

726: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
727: {
729:   lim->ops->view    = PetscLimiterView_VanLeer;
730:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
731:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
732:   return(0);
733: }

735: /*MC
736:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

738:   Level: intermediate

740: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
741: M*/

745: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
746: {
747:   PetscLimiter_VanLeer *l;
748:   PetscErrorCode    ierr;

752:   PetscNewLog(lim, &l);
753:   lim->data = l;

755:   PetscLimiterInitialize_VanLeer(lim);
756:   return(0);
757: }

761: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
762: {
763:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
764:   PetscErrorCode    ierr;

767:   PetscFree(l);
768:   return(0);
769: }

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

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

786: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
787: {
788:   PetscBool      iascii;

794:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
795:   if (iascii) {PetscLimiterView_VanAlbada_Ascii(lim, viewer);}
796:   return(0);
797: }

801: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
802: {
804:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
805:   return(0);
806: }

810: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
811: {
813:   lim->ops->view    = PetscLimiterView_VanAlbada;
814:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
815:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
816:   return(0);
817: }

819: /*MC
820:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

822:   Level: intermediate

824: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
825: M*/

829: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
830: {
831:   PetscLimiter_VanAlbada *l;
832:   PetscErrorCode    ierr;

836:   PetscNewLog(lim, &l);
837:   lim->data = l;

839:   PetscLimiterInitialize_VanAlbada(lim);
840:   return(0);
841: }

845: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
846: {
847:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
848:   PetscErrorCode    ierr;

851:   PetscFree(l);
852:   return(0);
853: }

857: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
858: {
859:   PetscViewerFormat format;
860:   PetscErrorCode    ierr;

863:   PetscViewerGetFormat(viewer, &format);
864:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
865:   return(0);
866: }

870: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
871: {
872:   PetscBool      iascii;

878:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
879:   if (iascii) {PetscLimiterView_Superbee_Ascii(lim, viewer);}
880:   return(0);
881: }

885: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
886: {
888:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
889:   return(0);
890: }

894: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
895: {
897:   lim->ops->view    = PetscLimiterView_Superbee;
898:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
899:   lim->ops->limit   = PetscLimiterLimit_Superbee;
900:   return(0);
901: }

903: /*MC
904:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

906:   Level: intermediate

908: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
909: M*/

913: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
914: {
915:   PetscLimiter_Superbee *l;
916:   PetscErrorCode    ierr;

920:   PetscNewLog(lim, &l);
921:   lim->data = l;

923:   PetscLimiterInitialize_Superbee(lim);
924:   return(0);
925: }

929: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
930: {
931:   PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
932:   PetscErrorCode    ierr;

935:   PetscFree(l);
936:   return(0);
937: }

941: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
942: {
943:   PetscViewerFormat format;
944:   PetscErrorCode    ierr;

947:   PetscViewerGetFormat(viewer, &format);
948:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
949:   return(0);
950: }

954: PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
955: {
956:   PetscBool      iascii;

962:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
963:   if (iascii) {PetscLimiterView_MC_Ascii(lim, viewer);}
964:   return(0);
965: }

969: /* aka Barth-Jespersen */
970: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
971: {
973:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
974:   return(0);
975: }

979: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
980: {
982:   lim->ops->view    = PetscLimiterView_MC;
983:   lim->ops->destroy = PetscLimiterDestroy_MC;
984:   lim->ops->limit   = PetscLimiterLimit_MC;
985:   return(0);
986: }

988: /*MC
989:   PETSCLIMITERMC = "mc" - A PetscLimiter object

991:   Level: intermediate

993: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
994: M*/

998: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
999: {
1000:   PetscLimiter_MC *l;
1001:   PetscErrorCode    ierr;

1005:   PetscNewLog(lim, &l);
1006:   lim->data = l;

1008:   PetscLimiterInitialize_MC(lim);
1009:   return(0);
1010: }

1012: PetscClassId PETSCFV_CLASSID = 0;

1014: PetscFunctionList PetscFVList              = NULL;
1015: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

1019: /*@C
1020:   PetscFVRegister - Adds a new PetscFV implementation

1022:   Not Collective

1024:   Input Parameters:
1025: + name        - The name of a new user-defined creation routine
1026: - create_func - The creation routine itself

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

1031:   Sample usage:
1032: .vb
1033:     PetscFVRegister("my_fv", MyPetscFVCreate);
1034: .ve

1036:   Then, your PetscFV type can be chosen with the procedural interface via
1037: .vb
1038:     PetscFVCreate(MPI_Comm, PetscFV *);
1039:     PetscFVSetType(PetscFV, "my_fv");
1040: .ve
1041:    or at runtime via the option
1042: .vb
1043:     -petscfv_type my_fv
1044: .ve

1046:   Level: advanced

1048: .keywords: PetscFV, register
1049: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()

1051: @*/
1052: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
1053: {

1057:   PetscFunctionListAdd(&PetscFVList, sname, function);
1058:   return(0);
1059: }

1063: /*@C
1064:   PetscFVSetType - Builds a particular PetscFV

1066:   Collective on PetscFV

1068:   Input Parameters:
1069: + fvm  - The PetscFV object
1070: - name - The kind of FVM space

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

1075:   Level: intermediate

1077: .keywords: PetscFV, set, type
1078: .seealso: PetscFVGetType(), PetscFVCreate()
1079: @*/
1080: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
1081: {
1082:   PetscErrorCode (*r)(PetscFV);
1083:   PetscBool      match;

1088:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1089:   if (match) return(0);

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

1095:   if (fvm->ops->destroy) {
1096:     (*fvm->ops->destroy)(fvm);
1097:     fvm->ops->destroy = NULL;
1098:   }
1099:   (*r)(fvm);
1100:   PetscObjectChangeTypeName((PetscObject) fvm, name);
1101:   return(0);
1102: }

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

1109:   Not Collective

1111:   Input Parameter:
1112: . fvm  - The PetscFV

1114:   Output Parameter:
1115: . name - The PetscFV type name

1117:   Level: intermediate

1119: .keywords: PetscFV, get, type, name
1120: .seealso: PetscFVSetType(), PetscFVCreate()
1121: @*/
1122: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1123: {

1129:   PetscFVRegisterAll();
1130:   *name = ((PetscObject) fvm)->type_name;
1131:   return(0);
1132: }

1136: /*@C
1137:   PetscFVView - Views a PetscFV

1139:   Collective on PetscFV

1141:   Input Parameter:
1142: + fvm - the PetscFV object to view
1143: - v   - the viewer

1145:   Level: developer

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

1155:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1156:   if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1157:   return(0);
1158: }

1162: /*@
1163:   PetscFVSetFromOptions - sets parameters in a PetscFV from the options database

1165:   Collective on PetscFV

1167:   Input Parameter:
1168: . fvm - the PetscFV object to set options for

1170:   Level: developer

1172: .seealso: PetscFVView()
1173: @*/
1174: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1175: {
1176:   const char    *defaultType;
1177:   char           name[256];
1178:   PetscBool      flg;

1183:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1184:   else                                 defaultType = ((PetscObject) fvm)->type_name;
1185:   PetscFVRegisterAll();

1187:   PetscObjectOptionsBegin((PetscObject) fvm);
1188:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1189:   if (flg) {
1190:     PetscFVSetType(fvm, name);
1191:   } else if (!((PetscObject) fvm)->type_name) {
1192:     PetscFVSetType(fvm, defaultType);
1193:   }
1194:   if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1195:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1196:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1197:   PetscLimiterSetFromOptions(fvm->limiter);
1198:   PetscOptionsEnd();
1199:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1200:   return(0);
1201: }

1205: /*@
1206:   PetscFVSetUp - Construct data structures for the PetscFV

1208:   Collective on PetscFV

1210:   Input Parameter:
1211: . fvm - the PetscFV object to setup

1213:   Level: developer

1215: .seealso: PetscFVView(), PetscFVDestroy()
1216: @*/
1217: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1218: {

1223:   PetscLimiterSetUp(fvm->limiter);
1224:   if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1225:   return(0);
1226: }

1230: /*@
1231:   PetscFVDestroy - Destroys a PetscFV object

1233:   Collective on PetscFV

1235:   Input Parameter:
1236: . fvm - the PetscFV object to destroy

1238:   Level: developer

1240: .seealso: PetscFVView()
1241: @*/
1242: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1243: {
1244:   PetscInt       i;

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

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

1254:   for (i = 0; i < (*fvm)->numComponents; i++) {
1255:     PetscFree((*fvm)->componentNames[i]);
1256:   }
1257:   PetscFree((*fvm)->componentNames);
1258:   PetscLimiterDestroy(&(*fvm)->limiter);
1259:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1260:   PetscFree((*fvm)->fluxWork);
1261:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1262:   PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);

1264:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1265:   PetscHeaderDestroy(fvm);
1266:   return(0);
1267: }

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

1274:   Collective on MPI_Comm

1276:   Input Parameter:
1277: . comm - The communicator for the PetscFV object

1279:   Output Parameter:
1280: . fvm - The PetscFV object

1282:   Level: beginner

1284: .seealso: PetscFVSetType(), PETSCFVUPWIND
1285: @*/
1286: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1287: {
1288:   PetscFV        f;

1293:   *fvm = NULL;
1294:   PetscFVInitializePackage();

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

1299:   PetscLimiterCreate(comm, &f->limiter);
1300:   f->numComponents    = 1;
1301:   f->dim              = 0;
1302:   f->computeGradients = PETSC_FALSE;
1303:   f->fluxWork         = NULL;
1304:   PetscCalloc1(f->numComponents,&f->componentNames);

1306:   *fvm = f;
1307:   return(0);
1308: }

1312: /*@
1313:   PetscFVSetLimiter - Set the limiter object

1315:   Logically collective on PetscFV

1317:   Input Parameters:
1318: + fvm - the PetscFV object
1319: - lim - The PetscLimiter

1321:   Level: developer

1323: .seealso: PetscFVGetLimiter()
1324: @*/
1325: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1326: {

1332:   PetscLimiterDestroy(&fvm->limiter);
1333:   PetscObjectReference((PetscObject) lim);
1334:   fvm->limiter = lim;
1335:   return(0);
1336: }

1340: /*@
1341:   PetscFVGetLimiter - Get the limiter object

1343:   Not collective

1345:   Input Parameter:
1346: . fvm - the PetscFV object

1348:   Output Parameter:
1349: . lim - The PetscLimiter

1351:   Level: developer

1353: .seealso: PetscFVSetLimiter()
1354: @*/
1355: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1356: {
1360:   *lim = fvm->limiter;
1361:   return(0);
1362: }

1366: /*@
1367:   PetscFVSetNumComponents - Set the number of field components

1369:   Logically collective on PetscFV

1371:   Input Parameters:
1372: + fvm - the PetscFV object
1373: - comp - The number of components

1375:   Level: developer

1377: .seealso: PetscFVGetNumComponents()
1378: @*/
1379: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1380: {

1385:   if (fvm->numComponents != comp) {
1386:     PetscInt i;

1388:     for (i = 0; i < fvm->numComponents; i++) {
1389:       PetscFree(fvm->componentNames[i]);
1390:     }
1391:     PetscFree(fvm->componentNames);
1392:     PetscCalloc1(comp,&fvm->componentNames);
1393:   }
1394:   fvm->numComponents = comp;
1395:   PetscFree(fvm->fluxWork);
1396:   PetscMalloc1(comp, &fvm->fluxWork);
1397:   return(0);
1398: }

1402: /*@
1403:   PetscFVGetNumComponents - Get the number of field components

1405:   Not collective

1407:   Input Parameter:
1408: . fvm - the PetscFV object

1410:   Output Parameter:
1411: , comp - The number of components

1413:   Level: developer

1415: .seealso: PetscFVSetNumComponents()
1416: @*/
1417: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1418: {
1422:   *comp = fvm->numComponents;
1423:   return(0);
1424: }

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

1431:   Logically collective on PetscFV
1432:   Input Parameters:
1433: + fvm - the PetscFV object
1434: . comp - the component number
1435: - name - the component name

1437:   Level: developer

1439: .seealso: PetscFVGetComponentName()
1440: @*/
1441: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1442: {

1446:   PetscFree(fvm->componentNames[comp]);
1447:   PetscStrallocpy(name,&fvm->componentNames[comp]);
1448:   return(0);
1449: }

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

1456:   Logically collective on PetscFV
1457:   Input Parameters:
1458: + fvm - the PetscFV object
1459: - comp - the component number

1461:   Output Parameter:
1462: . name - the component name

1464:   Level: developer

1466: .seealso: PetscFVSetComponentName()
1467: @*/
1468: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1469: {
1471:   *name = fvm->componentNames[comp];
1472:   return(0);
1473: }

1477: /*@
1478:   PetscFVSetSpatialDimension - Set the spatial dimension

1480:   Logically collective on PetscFV

1482:   Input Parameters:
1483: + fvm - the PetscFV object
1484: - dim - The spatial dimension

1486:   Level: developer

1488: .seealso: PetscFVGetSpatialDimension()
1489: @*/
1490: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1491: {
1494:   fvm->dim = dim;
1495:   return(0);
1496: }

1500: /*@
1501:   PetscFVGetSpatialDimension - Get the spatial dimension

1503:   Logically collective on PetscFV

1505:   Input Parameter:
1506: . fvm - the PetscFV object

1508:   Output Parameter:
1509: . dim - The spatial dimension

1511:   Level: developer

1513: .seealso: PetscFVSetSpatialDimension()
1514: @*/
1515: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1516: {
1520:   *dim = fvm->dim;
1521:   return(0);
1522: }

1526: /*@
1527:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1529:   Logically collective on PetscFV

1531:   Input Parameters:
1532: + fvm - the PetscFV object
1533: - computeGradients - Flag to compute cell gradients

1535:   Level: developer

1537: .seealso: PetscFVGetComputeGradients()
1538: @*/
1539: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1540: {
1543:   fvm->computeGradients = computeGradients;
1544:   return(0);
1545: }

1549: /*@
1550:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1552:   Not collective

1554:   Input Parameter:
1555: . fvm - the PetscFV object

1557:   Output Parameter:
1558: . computeGradients - Flag to compute cell gradients

1560:   Level: developer

1562: .seealso: PetscFVSetComputeGradients()
1563: @*/
1564: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1565: {
1569:   *computeGradients = fvm->computeGradients;
1570:   return(0);
1571: }

1575: /*@
1576:   PetscFVSetQuadrature - Set the quadrature object

1578:   Logically collective on PetscFV

1580:   Input Parameters:
1581: + fvm - the PetscFV object
1582: - q - The PetscQuadrature

1584:   Level: developer

1586: .seealso: PetscFVGetQuadrature()
1587: @*/
1588: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1589: {

1594:   PetscQuadratureDestroy(&fvm->quadrature);
1595:   PetscObjectReference((PetscObject) q);
1596:   fvm->quadrature = q;
1597:   return(0);
1598: }

1602: /*@
1603:   PetscFVGetQuadrature - Get the quadrature object

1605:   Not collective

1607:   Input Parameter:
1608: . fvm - the PetscFV object

1610:   Output Parameter:
1611: . lim - The PetscQuadrature

1613:   Level: developer

1615: .seealso: PetscFVSetQuadrature()
1616: @*/
1617: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1618: {
1622:   if (!fvm->quadrature) {
1623:     /* Create default 1-point quadrature */
1624:     PetscReal     *points, *weights;

1627:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1628:     PetscCalloc1(fvm->dim, &points);
1629:     PetscMalloc1(1, &weights);
1630:     weights[0] = 1.0;
1631:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, points, weights);
1632:   }
1633:   *q = fvm->quadrature;
1634:   return(0);
1635: }

1639: /*@
1640:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1642:   Not collective

1644:   Input Parameter:
1645: . fvm - The PetscFV object

1647:   Output Parameter:
1648: . sp - The PetscDualSpace object

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

1652:   Level: developer

1654: .seealso: PetscFVCreate()
1655: @*/
1656: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1657: {
1661:   if (!fvm->dualSpace) {
1662:     DM              K;
1663:     PetscQuadrature q;
1664:     PetscInt        dim;
1665:     PetscErrorCode  ierr;

1667:     PetscFVGetSpatialDimension(fvm, &dim);
1668:     PetscFVGetQuadrature(fvm, &q);
1669:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1670:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1671:     PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1672:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1673:     DMDestroy(&K);
1674:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, 1);
1675:     PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, 0, q);
1676:   }
1677:   *sp = fvm->dualSpace;
1678:   return(0);
1679: }

1683: /*@
1684:   PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product

1686:   Not collective

1688:   Input Parameters:
1689: + fvm - The PetscFV object
1690: - sp  - The PetscDualSpace object

1692:   Level: intermediate

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

1696: .seealso: PetscFVCreate()
1697: @*/
1698: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1699: {

1705:   PetscDualSpaceDestroy(&fvm->dualSpace);
1706:   fvm->dualSpace = sp;
1707:   PetscObjectReference((PetscObject) fvm->dualSpace);
1708:   return(0);
1709: }

1713: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1714: {
1715:   PetscInt         npoints;
1716:   const PetscReal *points;
1717:   PetscErrorCode   ierr;

1724:   PetscQuadratureGetData(fvm->quadrature, NULL, &npoints, &points, NULL);
1725:   if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1726:   if (B) *B = fvm->B;
1727:   if (D) *D = fvm->D;
1728:   if (H) *H = fvm->H;
1729:   return(0);
1730: }

1734: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1735: {
1736:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1737:   PetscInt         dim;      /* Spatial dimension */
1738:   PetscInt         comp;     /* Field components */
1739:   PetscInt         p, d, c, e;
1740:   PetscErrorCode   ierr;

1748:   PetscFVGetSpatialDimension(fvm, &dim);
1749:   PetscFVGetNumComponents(fvm, &comp);
1750:   if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1751:   if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1752:   if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1753:   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;}
1754:   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;}
1755:   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;}
1756:   return(0);
1757: }

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

1767:   if (B && *B) {PetscFree(*B);}
1768:   if (D && *D) {PetscFree(*D);}
1769:   if (H && *H) {PetscFree(*H);}
1770:   return(0);
1771: }

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

1778:   Input Parameters:
1779: + fvm      - The PetscFV object
1780: . numFaces - The number of cell faces which are not constrained
1781: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1783:   Level: developer

1785: .seealso: PetscFVCreate()
1786: @*/
1787: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1788: {

1793:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1794:   return(0);
1795: }

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

1802:   Not collective

1804:   Input Parameters:
1805: + fvm          - The PetscFV object for the field being integrated
1806: . prob         - The PetscDS specifing the discretizations and continuum functions
1807: . field        - The field being integrated
1808: . Nf           - The number of faces in the chunk
1809: . fgeom        - The face geometry for each face in the chunk
1810: . neighborVol  - The volume for each pair of cells in the chunk
1811: . uL           - The state from the cell on the left
1812: - uR           - The state from the cell on the right

1814:   Output Parameter
1815: + fluxL        - the left fluxes for each face
1816: - fluxR        - the right fluxes for each face
1817: */
1818: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1819:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1820: {

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

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

1836:   Input Parameter:
1837: . fv - The initial PetscFV

1839:   Output Parameter:
1840: . fvRef - The refined PetscFV

1842:   Level: developer

1844: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1845: @*/
1846: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1847: {
1848:   PetscDualSpace   Q, Qref;
1849:   DM               K, Kref;
1850:   PetscQuadrature  q, qref;
1851:   CellRefiner      cellRefiner;
1852:   PetscReal       *v0;
1853:   PetscReal       *jac, *invjac;
1854:   PetscInt         numComp, numSubelements, s;
1855:   PetscErrorCode   ierr;

1858:   PetscFVGetDualSpace(fv, &Q);
1859:   PetscFVGetQuadrature(fv, &q);
1860:   PetscDualSpaceGetDM(Q, &K);
1861:   /* Create dual space */
1862:   PetscDualSpaceDuplicate(Q, &Qref);
1863:   DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1864:   PetscDualSpaceSetDM(Qref, Kref);
1865:   DMDestroy(&Kref);
1866:   PetscDualSpaceSetUp(Qref);
1867:   /* Create volume */
1868:   PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1869:   PetscFVSetDualSpace(*fvRef, Qref);
1870:   PetscFVGetNumComponents(fv,    &numComp);
1871:   PetscFVSetNumComponents(*fvRef, numComp);
1872:   PetscFVSetUp(*fvRef);
1873:   /* Create quadrature */
1874:   DMPlexGetCellRefiner_Internal(K, &cellRefiner);
1875:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1876:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1877:   PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1878:   for (s = 0; s < numSubelements; ++s) {
1879:     PetscQuadrature  qs;
1880:     const PetscReal *points, *weights;
1881:     PetscReal       *p, *w;
1882:     PetscInt         dim, npoints, np;

1884:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1885:     PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
1886:     np   = npoints/numSubelements;
1887:     PetscMalloc1(np*dim,&p);
1888:     PetscMalloc1(np,&w);
1889:     PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1890:     PetscMemcpy(w, &weights[s*np],    np     * sizeof(PetscReal));
1891:     PetscQuadratureSetData(qs, dim, np, p, w);
1892:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1893:     PetscQuadratureDestroy(&qs);
1894:   }
1895:   CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1896:   PetscFVSetQuadrature(*fvRef, qref);
1897:   PetscQuadratureDestroy(&qref);
1898:   PetscDualSpaceDestroy(&Qref);
1899:   return(0);
1900: }

1904: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1905: {
1906:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1910:   PetscFree(b);
1911:   return(0);
1912: }

1916: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1917: {
1918:   PetscInt          Nc, c;
1919:   PetscViewerFormat format;
1920:   PetscErrorCode    ierr;

1923:   PetscFVGetNumComponents(fv, &Nc);
1924:   PetscViewerGetFormat(viewer, &format);
1925:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1926:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1927:   for (c = 0; c < Nc; c++) {
1928:     if (fv->componentNames[c]) {
1929:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1930:     }
1931:   }
1932:   return(0);
1933: }

1937: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1938: {
1939:   PetscBool      iascii;

1945:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1946:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1947:   return(0);
1948: }

1952: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1953: {
1955:   return(0);
1956: }

1960: /*
1961:   neighborVol[f*2+0] contains the left  geom
1962:   neighborVol[f*2+1] contains the right geom
1963: */
1964: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1965:                                                   PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1966: {
1967:   void         (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
1968:   void          *rctx;
1969:   PetscScalar   *flux = fvm->fluxWork;
1970:   PetscInt       dim, pdim, totDim, Nc, off, f, d;

1974:   PetscDSGetTotalComponents(prob, &Nc);
1975:   PetscDSGetTotalDimension(prob, &totDim);
1976:   PetscDSGetFieldOffset(prob, field, &off);
1977:   PetscDSGetRiemannSolver(prob, field, &riemann);
1978:   PetscDSGetContext(prob, field, &rctx);
1979:   PetscFVGetSpatialDimension(fvm, &dim);
1980:   PetscFVGetNumComponents(fvm, &pdim);
1981:   for (f = 0; f < Nf; ++f) {
1982:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
1983:     for (d = 0; d < pdim; ++d) {
1984:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1985:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1986:     }
1987:   }
1988:   return(0);
1989: }

1993: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1994: {
1996:   fvm->ops->setfromoptions          = NULL;
1997:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1998:   fvm->ops->view                    = PetscFVView_Upwind;
1999:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
2000:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
2001:   return(0);
2002: }

2004: /*MC
2005:   PETSCFVUPWIND = "upwind" - A PetscFV object

2007:   Level: intermediate

2009: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2010: M*/

2014: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
2015: {
2016:   PetscFV_Upwind *b;
2017:   PetscErrorCode  ierr;

2021:   PetscNewLog(fvm,&b);
2022:   fvm->data = b;

2024:   PetscFVInitialize_Upwind(fvm);
2025:   return(0);
2026: }

2028: #include <petscblaslapack.h>

2032: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
2033: {
2034:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2035:   PetscErrorCode        ierr;

2038:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
2039:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2040:   PetscFree(ls);
2041:   return(0);
2042: }

2046: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
2047: {
2048:   PetscInt          Nc, c;
2049:   PetscViewerFormat format;
2050:   PetscErrorCode    ierr;

2053:   PetscFVGetNumComponents(fv, &Nc);
2054:   PetscViewerGetFormat(viewer, &format);
2055:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
2056:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
2057:   for (c = 0; c < Nc; c++) {
2058:     if (fv->componentNames[c]) {
2059:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
2060:     }
2061:   }
2062:   return(0);
2063: }

2067: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2068: {
2069:   PetscBool      iascii;

2075:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2076:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2077:   return(0);
2078: }

2082: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2083: {
2085:   return(0);
2086: }

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

2099:   if (debug) {
2100:     PetscMalloc1(m*n,&Aback);
2101:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2102:   }

2104:   PetscBLASIntCast(m,&M);
2105:   PetscBLASIntCast(n,&N);
2106:   PetscBLASIntCast(mstride,&lda);
2107:   PetscBLASIntCast(worksize,&ldwork);
2108:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2109:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2110:   PetscFPTrapPop();
2111:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2112:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2114:   /* Extract an explicit representation of Q */
2115:   Q    = Ainv;
2116:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
2117:   K    = N;                     /* full rank */
2118:   LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
2119:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

2121:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2122:   Alpha = 1.0;
2123:   ldb   = lda;
2124:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2125:   /* Ainv is Q, overwritten with inverse */

2127:   if (debug) {                      /* Check that pseudo-inverse worked */
2128:     PetscScalar  Beta = 0.0;
2129:     PetscBLASInt ldc;
2130:     K   = N;
2131:     ldc = N;
2132:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2133:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2134:     PetscFree(Aback);
2135:   }
2136:   return(0);
2137: }

2141: /* Overwrites A. Can handle degenerate problems and m<n. */
2142: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2143: {
2144:   PetscBool      debug = PETSC_FALSE;
2145:   PetscScalar   *Brhs, *Aback;
2146:   PetscScalar   *tmpwork;
2147:   PetscReal      rcond;
2148: #if defined (PETSC_USE_COMPLEX)
2149:   PetscInt       rworkSize;
2150:   PetscReal     *rwork;
2151: #endif
2152:   PetscInt       i, j, maxmn;
2153:   PetscBLASInt   M, N, lda, ldb, ldwork;
2154:   PetscBLASInt   nrhs, irank, info;

2158:   if (debug) {
2159:     PetscMalloc1(m*n,&Aback);
2160:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2161:   }

2163:   /* initialize to identity */
2164:   tmpwork = Ainv;
2165:   Brhs = work;
2166:   maxmn = PetscMax(m,n);
2167:   for (j=0; j<maxmn; j++) {
2168:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2169:   }

2171:   PetscBLASIntCast(m,&M);
2172:   PetscBLASIntCast(n,&N);
2173:   PetscBLASIntCast(mstride,&lda);
2174:   PetscBLASIntCast(maxmn,&ldb);
2175:   PetscBLASIntCast(worksize,&ldwork);
2176:   rcond = -1;
2177:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2178:   nrhs  = M;
2179: #if defined(PETSC_USE_COMPLEX)
2180:   rworkSize = 5 * PetscMin(M,N);
2181:   PetscMalloc1(rworkSize,&rwork);
2182:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2183:   PetscFPTrapPop();
2184:   PetscFree(rwork);
2185: #else
2186:   nrhs  = M;
2187:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2188:   PetscFPTrapPop();
2189: #endif
2190:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2191:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2192:   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");

2194:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2195:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2196:   for (i=0; i<n; i++) {
2197:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2198:   }
2199:   return(0);
2200: }

2202: #if 0
2205: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2206: {
2207:   PetscReal       grad[2] = {0, 0};
2208:   const PetscInt *faces;
2209:   PetscInt        numFaces, f;
2210:   PetscErrorCode  ierr;

2213:   DMPlexGetConeSize(dm, cell, &numFaces);
2214:   DMPlexGetCone(dm, cell, &faces);
2215:   for (f = 0; f < numFaces; ++f) {
2216:     const PetscInt *fcells;
2217:     const CellGeom *cg1;
2218:     const FaceGeom *fg;

2220:     DMPlexGetSupport(dm, faces[f], &fcells);
2221:     DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2222:     for (i = 0; i < 2; ++i) {
2223:       PetscScalar du;

2225:       if (fcells[i] == c) continue;
2226:       DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2227:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2228:       grad[0] += fg->grad[!i][0] * du;
2229:       grad[1] += fg->grad[!i][1] * du;
2230:     }
2231:   }
2232:   PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2233:   return(0);
2234: }
2235: #endif

2239: /*
2240:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

2242:   Input Parameters:
2243: + fvm      - The PetscFV object
2244: . numFaces - The number of cell faces which are not constrained
2245: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

2247:   Level: developer

2249: .seealso: PetscFVCreate()
2250: */
2251: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2252: {
2253:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *) fvm->data;
2254:   const PetscBool       useSVD   = PETSC_TRUE;
2255:   const PetscInt        maxFaces = ls->maxFaces;
2256:   PetscInt              dim, f, d;
2257:   PetscErrorCode        ierr;

2260:   if (numFaces > maxFaces) {
2261:     if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2262:     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2263:   }
2264:   PetscFVGetSpatialDimension(fvm, &dim);
2265:   for (f = 0; f < numFaces; ++f) {
2266:     for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2267:   }
2268:   /* Overwrites B with garbage, returns Binv in row-major format */
2269:   if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2270:   else        {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2271:   for (f = 0; f < numFaces; ++f) {
2272:     for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2273:   }
2274:   return(0);
2275: }

2279: /*
2280:   neighborVol[f*2+0] contains the left  geom
2281:   neighborVol[f*2+1] contains the right geom
2282: */
2283: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2284:                                                         PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2285: {
2286:   void         (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
2287:   void          *rctx;
2288:   PetscScalar   *flux = fvm->fluxWork;
2289:   PetscInt       dim, pdim, Nc, totDim, off, f, d;

2293:   PetscDSGetTotalComponents(prob, &Nc);
2294:   PetscDSGetTotalDimension(prob, &totDim);
2295:   PetscDSGetFieldOffset(prob, field, &off);
2296:   PetscDSGetRiemannSolver(prob, field, &riemann);
2297:   PetscDSGetContext(prob, field, &rctx);
2298:   PetscFVGetSpatialDimension(fvm, &dim);
2299:   PetscFVGetNumComponents(fvm, &pdim);
2300:   for (f = 0; f < Nf; ++f) {
2301:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
2302:     for (d = 0; d < pdim; ++d) {
2303:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2304:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2305:     }
2306:   }
2307:   return(0);
2308: }

2312: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2313: {
2314:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2315:   PetscInt              dim, m, n, nrhs, minwork;
2316:   PetscErrorCode        ierr;

2320:   PetscFVGetSpatialDimension(fvm, &dim);
2321:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2322:   ls->maxFaces = maxFaces;
2323:   m       = ls->maxFaces;
2324:   n       = dim;
2325:   nrhs    = ls->maxFaces;
2326:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2327:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2328:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2329:   return(0);
2330: }

2334: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2335: {
2337:   fvm->ops->setfromoptions          = NULL;
2338:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2339:   fvm->ops->view                    = PetscFVView_LeastSquares;
2340:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2341:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2342:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2343:   return(0);
2344: }

2346: /*MC
2347:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2349:   Level: intermediate

2351: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2352: M*/

2356: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2357: {
2358:   PetscFV_LeastSquares *ls;
2359:   PetscErrorCode        ierr;

2363:   PetscNewLog(fvm, &ls);
2364:   fvm->data = ls;

2366:   ls->maxFaces = -1;
2367:   ls->workSize = -1;
2368:   ls->B        = NULL;
2369:   ls->Binv     = NULL;
2370:   ls->tau      = NULL;
2371:   ls->work     = NULL;

2373:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2374:   PetscFVInitialize_LeastSquares(fvm);
2375:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2376:   return(0);
2377: }

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

2384:   Not collective

2386:   Input parameters:
2387: + fvm      - The PetscFV object
2388: - maxFaces - The maximum number of cell faces

2390:   Level: intermediate

2392: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2393: @*/
2394: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2395: {

2400:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2401:   return(0);
2402: }