Actual source code: dtfv.c

petsc-3.6.1 2015-07-22
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((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((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: {

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

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

1253:   PetscLimiterDestroy(&(*fvm)->limiter);
1254:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1255:   PetscFree((*fvm)->fluxWork);
1256:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1257:   PetscFVRestoreTabulation((*fvm), 0, NULL, &(*fvm)->B, &(*fvm)->D, NULL /*&(*fvm)->H*/);

1259:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1260:   PetscHeaderDestroy(fvm);
1261:   return(0);
1262: }

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

1269:   Collective on MPI_Comm

1271:   Input Parameter:
1272: . comm - The communicator for the PetscFV object

1274:   Output Parameter:
1275: . fvm - The PetscFV object

1277:   Level: beginner

1279: .seealso: PetscFVSetType(), PETSCFVUPWIND
1280: @*/
1281: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1282: {
1283:   PetscFV        f;

1288:   *fvm = NULL;
1289:   PetscFVInitializePackage();

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

1294:   PetscLimiterCreate(comm, &f->limiter);
1295:   f->numComponents    = 1;
1296:   f->dim              = 0;
1297:   f->computeGradients = PETSC_FALSE;
1298:   f->fluxWork         = NULL;

1300:   *fvm = f;
1301:   return(0);
1302: }

1306: /*@
1307:   PetscFVSetLimiter - Set the limiter object

1309:   Logically collective on PetscFV

1311:   Input Parameters:
1312: + fvm - the PetscFV object
1313: - lim - The PetscLimiter

1315:   Level: developer

1317: .seealso: PetscFVGetLimiter()
1318: @*/
1319: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1320: {

1326:   PetscLimiterDestroy(&fvm->limiter);
1327:   PetscObjectReference((PetscObject) lim);
1328:   fvm->limiter = lim;
1329:   return(0);
1330: }

1334: /*@
1335:   PetscFVGetLimiter - Get the limiter object

1337:   Not collective

1339:   Input Parameter:
1340: . fvm - the PetscFV object

1342:   Output Parameter:
1343: . lim - The PetscLimiter

1345:   Level: developer

1347: .seealso: PetscFVSetLimiter()
1348: @*/
1349: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1350: {
1354:   *lim = fvm->limiter;
1355:   return(0);
1356: }

1360: /*@
1361:   PetscFVSetNumComponents - Set the number of field components

1363:   Logically collective on PetscFV

1365:   Input Parameters:
1366: + fvm - the PetscFV object
1367: - comp - The number of components

1369:   Level: developer

1371: .seealso: PetscFVGetNumComponents()
1372: @*/
1373: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1374: {

1379:   fvm->numComponents = comp;
1380:   PetscFree(fvm->fluxWork);
1381:   PetscMalloc1(comp, &fvm->fluxWork);
1382:   return(0);
1383: }

1387: /*@
1388:   PetscFVGetNumComponents - Get the number of field components

1390:   Not collective

1392:   Input Parameter:
1393: . fvm - the PetscFV object

1395:   Output Parameter:
1396: , comp - The number of components

1398:   Level: developer

1400: .seealso: PetscFVSetNumComponents()
1401: @*/
1402: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1403: {
1407:   *comp = fvm->numComponents;
1408:   return(0);
1409: }

1413: /*@
1414:   PetscFVSetSpatialDimension - Set the spatial dimension

1416:   Logically collective on PetscFV

1418:   Input Parameters:
1419: + fvm - the PetscFV object
1420: - dim - The spatial dimension

1422:   Level: developer

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

1436: /*@
1437:   PetscFVGetSpatialDimension - Get the spatial dimension

1439:   Logically collective on PetscFV

1441:   Input Parameter:
1442: . fvm - the PetscFV object

1444:   Output Parameter:
1445: . dim - The spatial dimension

1447:   Level: developer

1449: .seealso: PetscFVSetSpatialDimension()
1450: @*/
1451: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1452: {
1456:   *dim = fvm->dim;
1457:   return(0);
1458: }

1462: /*@
1463:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1465:   Logically collective on PetscFV

1467:   Input Parameters:
1468: + fvm - the PetscFV object
1469: - computeGradients - Flag to compute cell gradients

1471:   Level: developer

1473: .seealso: PetscFVGetComputeGradients()
1474: @*/
1475: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1476: {
1479:   fvm->computeGradients = computeGradients;
1480:   return(0);
1481: }

1485: /*@
1486:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1488:   Not collective

1490:   Input Parameter:
1491: . fvm - the PetscFV object

1493:   Output Parameter:
1494: . computeGradients - Flag to compute cell gradients

1496:   Level: developer

1498: .seealso: PetscFVSetComputeGradients()
1499: @*/
1500: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1501: {
1505:   *computeGradients = fvm->computeGradients;
1506:   return(0);
1507: }

1511: /*@
1512:   PetscFVSetQuadrature - Set the quadrature object

1514:   Logically collective on PetscFV

1516:   Input Parameters:
1517: + fvm - the PetscFV object
1518: - q - The PetscQuadrature

1520:   Level: developer

1522: .seealso: PetscFVGetQuadrature()
1523: @*/
1524: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1525: {

1530:   PetscQuadratureDestroy(&fvm->quadrature);
1531:   PetscObjectReference((PetscObject) q);
1532:   fvm->quadrature = q;
1533:   return(0);
1534: }

1538: /*@
1539:   PetscFVGetQuadrature - Get the quadrature object

1541:   Not collective

1543:   Input Parameter:
1544: . fvm - the PetscFV object

1546:   Output Parameter:
1547: . lim - The PetscQuadrature

1549:   Level: developer

1551: .seealso: PetscFVSetQuadrature()
1552: @*/
1553: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1554: {
1558:   if (!fvm->quadrature) {
1559:     /* Create default 1-point quadrature */
1560:     PetscReal     *points, *weights;

1563:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1564:     PetscCalloc1(fvm->dim, &points);
1565:     PetscMalloc1(1, &weights);
1566:     weights[0] = 1.0;
1567:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, points, weights);
1568:   }
1569:   *q = fvm->quadrature;
1570:   return(0);
1571: }

1575: /*@
1576:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1578:   Not collective

1580:   Input Parameter:
1581: . fvm - The PetscFV object

1583:   Output Parameter:
1584: . sp - The PetscDualSpace object

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

1588:   Level: developer

1590: .seealso: PetscFVCreate()
1591: @*/
1592: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1593: {
1597:   if (!fvm->dualSpace) {
1598:     DM              K;
1599:     PetscQuadrature q;
1600:     PetscInt        dim;
1601:     PetscErrorCode  ierr;

1603:     PetscFVGetSpatialDimension(fvm, &dim);
1604:     PetscFVGetQuadrature(fvm, &q);
1605:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1606:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1607:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, 1);
1608:     PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, 0, q);
1609:     PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K); /* TODO: The reference cell type should be held by the discretization object */
1610:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1611:     DMDestroy(&K);
1612:   }
1613:   *sp = fvm->dualSpace;
1614:   return(0);
1615: }

1619: /*@
1620:   PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product

1622:   Not collective

1624:   Input Parameters:
1625: + fvm - The PetscFV object
1626: - sp  - The PetscDualSpace object

1628:   Level: intermediate

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

1632: .seealso: PetscFVCreate()
1633: @*/
1634: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1635: {

1641:   PetscDualSpaceDestroy(&fvm->dualSpace);
1642:   fvm->dualSpace = sp;
1643:   PetscObjectReference((PetscObject) fvm->dualSpace);
1644:   return(0);
1645: }

1649: PetscErrorCode PetscFVGetDefaultTabulation(PetscFV fvm, PetscReal **B, PetscReal **D, PetscReal **H)
1650: {
1651:   PetscInt         npoints;
1652:   const PetscReal *points;
1653:   PetscErrorCode   ierr;

1660:   PetscQuadratureGetData(fvm->quadrature, NULL, &npoints, &points, NULL);
1661:   if (!fvm->B) {PetscFVGetTabulation(fvm, npoints, points, &fvm->B, &fvm->D, NULL/*&fvm->H*/);}
1662:   if (B) *B = fvm->B;
1663:   if (D) *D = fvm->D;
1664:   if (H) *H = fvm->H;
1665:   return(0);
1666: }

1670: PetscErrorCode PetscFVGetTabulation(PetscFV fvm, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
1671: {
1672:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1673:   PetscInt         dim;      /* Spatial dimension */
1674:   PetscInt         comp;     /* Field components */
1675:   PetscInt         p, d, c, e;
1676:   PetscErrorCode   ierr;

1684:   PetscFVGetSpatialDimension(fvm, &dim);
1685:   PetscFVGetNumComponents(fvm, &comp);
1686:   if (B) {PetscMalloc1(npoints*pdim*comp, B);}
1687:   if (D) {PetscMalloc1(npoints*pdim*comp*dim, D);}
1688:   if (H) {PetscMalloc1(npoints*pdim*comp*dim*dim, H);}
1689:   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;}
1690:   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;}
1691:   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;}
1692:   return(0);
1693: }

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

1703:   if (B && *B) {PetscFree(*B);}
1704:   if (D && *D) {PetscFree(*D);}
1705:   if (H && *H) {PetscFree(*H);}
1706:   return(0);
1707: }

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

1714:   Input Parameters:
1715: + fvm      - The PetscFV object
1716: . numFaces - The number of cell faces which are not constrained
1717: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1719:   Level: developer

1721: .seealso: PetscFVCreate()
1722: @*/
1723: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1724: {

1729:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1730:   return(0);
1731: }

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

1738:   Not collective

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

1750:   Output Parameter
1751: + fluxL        - the left fluxes for each face
1752: - fluxR        - the right fluxes for each face
1753: */
1754: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1755:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1756: {

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

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

1772:   Input Parameter:
1773: . fv - The initial PetscFV

1775:   Output Parameter:
1776: . fvRef - The refined PetscFV

1778:   Level: developer

1780: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1781: @*/
1782: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1783: {
1784:   PetscDualSpace   Q, Qref;
1785:   DM               K, Kref;
1786:   PetscQuadrature  q, qref;
1787:   CellRefiner      cellRefiner;
1788:   PetscReal       *v0;
1789:   PetscReal       *jac, *invjac;
1790:   PetscInt         numComp, numSubelements, s;
1791:   PetscErrorCode   ierr;

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

1820:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1821:     PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
1822:     np   = npoints/numSubelements;
1823:     PetscMalloc1(np*dim,&p);
1824:     PetscMalloc1(np,&w);
1825:     PetscMemcpy(p, &points[s*np*dim], np*dim * sizeof(PetscReal));
1826:     PetscMemcpy(w, &weights[s*np],    np     * sizeof(PetscReal));
1827:     PetscQuadratureSetData(qs, dim, np, p, w);
1828:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1829:     PetscQuadratureDestroy(&qs);
1830:   }
1831:   CellRefinerRestoreAffineTransforms_Internal(cellRefiner, &numSubelements, &v0, &jac, &invjac);
1832:   PetscFVSetQuadrature(*fvRef, qref);
1833:   PetscQuadratureDestroy(&qref);
1834:   PetscDualSpaceDestroy(&Qref);
1835:   return(0);
1836: }

1840: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1841: {
1842:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1846:   PetscFree(b);
1847:   return(0);
1848: }

1852: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1853: {
1854:   PetscInt          Nc;
1855:   PetscViewerFormat format;
1856:   PetscErrorCode    ierr;

1859:   PetscFVGetNumComponents(fv, &Nc);
1860:   PetscViewerGetFormat(viewer, &format);
1861:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1862:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1863:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1864:   } else {
1865:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1866:   }
1867:   return(0);
1868: }

1872: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1873: {
1874:   PetscBool      iascii;

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

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

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

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

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

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

1942:   Level: intermediate

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

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

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

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

1963: #include <petscblaslapack.h>

1967: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1968: {
1969:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1970:   PetscErrorCode        ierr;

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

1981: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1982: {
1983:   PetscInt          Nc;
1984:   PetscViewerFormat format;
1985:   PetscErrorCode    ierr;

1988:   PetscFVGetNumComponents(fv, &Nc);
1989:   PetscViewerGetFormat(viewer, &format);
1990:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1991:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1992:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1993:   } else {
1994:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1995:   }
1996:   return(0);
1997: }

2001: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
2002: {
2003:   PetscBool      iascii;

2009:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2010:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
2011:   return(0);
2012: }

2016: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2017: {
2019:   return(0);
2020: }

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

2033:   if (debug) {
2034:     PetscMalloc1(m*n,&Aback);
2035:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2036:   }

2038:   PetscBLASIntCast(m,&M);
2039:   PetscBLASIntCast(n,&N);
2040:   PetscBLASIntCast(mstride,&lda);
2041:   PetscBLASIntCast(worksize,&ldwork);
2042:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2043:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2044:   PetscFPTrapPop();
2045:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2046:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

2048:   /* Extract an explicit representation of Q */
2049:   Q    = Ainv;
2050:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
2051:   K    = N;                     /* full rank */
2052:   LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
2053:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

2055:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2056:   Alpha = 1.0;
2057:   ldb   = lda;
2058:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2059:   /* Ainv is Q, overwritten with inverse */

2061:   if (debug) {                      /* Check that pseudo-inverse worked */
2062:     PetscScalar  Beta = 0.0;
2063:     PetscBLASInt ldc;
2064:     K   = N;
2065:     ldc = N;
2066:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2067:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
2068:     PetscFree(Aback);
2069:   }
2070:   return(0);
2071: }

2075: /* Overwrites A. Can handle degenerate problems and m<n. */
2076: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2077: {
2078:   PetscBool      debug = PETSC_FALSE;
2079:   PetscScalar   *Brhs, *Aback;
2080:   PetscScalar   *tmpwork;
2081:   PetscReal      rcond;
2082:   PetscInt       i, j, maxmn;
2083:   PetscBLASInt   M, N, lda, ldb, ldwork;
2084: #if !defined(PETSC_USE_COMPLEX)
2085:   PetscBLASInt   nrhs, irank, info;
2086: #endif

2090:   if (debug) {
2091:     PetscMalloc1(m*n,&Aback);
2092:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
2093:   }

2095:   /* initialize to identity */
2096:   tmpwork = Ainv;
2097:   Brhs = work;
2098:   maxmn = PetscMax(m,n);
2099:   for (j=0; j<maxmn; j++) {
2100:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2101:   }

2103:   PetscBLASIntCast(m,&M);
2104:   PetscBLASIntCast(n,&N);
2105:   PetscBLASIntCast(mstride,&lda);
2106:   PetscBLASIntCast(maxmn,&ldb);
2107:   PetscBLASIntCast(worksize,&ldwork);
2108:   rcond = -1;
2109:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
2110: #if defined(PETSC_USE_COMPLEX)
2111:   if (tmpwork && rcond) rcond = 0.0; /* Get rid of compiler warning */
2112:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "I don't think this makes sense for complex numbers");
2113: #else
2114:   nrhs  = M;
2115:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2116:   PetscFPTrapPop();
2117:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2118:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2119:   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");

2121:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
2122:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
2123:   for (i=0; i<n; i++) {
2124:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
2125:   }
2126:   return(0);
2127: #endif
2128: }

2130: #if 0
2133: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2134: {
2135:   PetscReal       grad[2] = {0, 0};
2136:   const PetscInt *faces;
2137:   PetscInt        numFaces, f;
2138:   PetscErrorCode  ierr;

2141:   DMPlexGetConeSize(dm, cell, &numFaces);
2142:   DMPlexGetCone(dm, cell, &faces);
2143:   for (f = 0; f < numFaces; ++f) {
2144:     const PetscInt *fcells;
2145:     const CellGeom *cg1;
2146:     const FaceGeom *fg;

2148:     DMPlexGetSupport(dm, faces[f], &fcells);
2149:     DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
2150:     for (i = 0; i < 2; ++i) {
2151:       PetscScalar du;

2153:       if (fcells[i] == c) continue;
2154:       DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
2155:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2156:       grad[0] += fg->grad[!i][0] * du;
2157:       grad[1] += fg->grad[!i][1] * du;
2158:     }
2159:   }
2160:   PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
2161:   return(0);
2162: }
2163: #endif

2167: /*
2168:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

2170:   Input Parameters:
2171: + fvm      - The PetscFV object
2172: . numFaces - The number of cell faces which are not constrained
2173: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

2175:   Level: developer

2177: .seealso: PetscFVCreate()
2178: */
2179: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2180: {
2181:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *) fvm->data;
2182:   const PetscBool       useSVD   = PETSC_TRUE;
2183:   const PetscInt        maxFaces = ls->maxFaces;
2184:   PetscInt              dim, f, d;
2185:   PetscErrorCode        ierr;

2188:   if (numFaces > maxFaces) {
2189:     if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2190:     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2191:   }
2192:   PetscFVGetSpatialDimension(fvm, &dim);
2193:   for (f = 0; f < numFaces; ++f) {
2194:     for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2195:   }
2196:   /* Overwrites B with garbage, returns Binv in row-major format */
2197:   if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2198:   else        {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
2199:   for (f = 0; f < numFaces; ++f) {
2200:     for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
2201:   }
2202:   return(0);
2203: }

2207: /*
2208:   neighborVol[f*2+0] contains the left  geom
2209:   neighborVol[f*2+1] contains the right geom
2210: */
2211: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2212:                                                         PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2213: {
2214:   void         (*riemann)(PetscInt, PetscInt, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx);
2215:   void          *rctx;
2216:   PetscScalar   *flux = fvm->fluxWork;
2217:   PetscInt       dim, pdim, Nc, totDim, off, f, d;

2221:   PetscDSGetTotalComponents(prob, &Nc);
2222:   PetscDSGetTotalDimension(prob, &totDim);
2223:   PetscDSGetFieldOffset(prob, field, &off);
2224:   PetscDSGetRiemannSolver(prob, field, &riemann);
2225:   PetscDSGetContext(prob, field, &rctx);
2226:   PetscFVGetSpatialDimension(fvm, &dim);
2227:   PetscFVGetNumComponents(fvm, &pdim);
2228:   for (f = 0; f < Nf; ++f) {
2229:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], flux, rctx);
2230:     for (d = 0; d < pdim; ++d) {
2231:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2232:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2233:     }
2234:   }
2235:   return(0);
2236: }

2240: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2241: {
2242:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2243:   PetscInt              dim, m, n, nrhs, minwork;
2244:   PetscErrorCode        ierr;

2248:   PetscFVGetSpatialDimension(fvm, &dim);
2249:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2250:   ls->maxFaces = maxFaces;
2251:   m       = ls->maxFaces;
2252:   n       = dim;
2253:   nrhs    = ls->maxFaces;
2254:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2255:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2256:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2257:   return(0);
2258: }

2262: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2263: {
2265:   fvm->ops->setfromoptions          = NULL;
2266:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2267:   fvm->ops->view                    = PetscFVView_LeastSquares;
2268:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2269:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2270:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2271:   return(0);
2272: }

2274: /*MC
2275:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2277:   Level: intermediate

2279: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2280: M*/

2284: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2285: {
2286:   PetscFV_LeastSquares *ls;
2287:   PetscErrorCode        ierr;

2291:   PetscNewLog(fvm, &ls);
2292:   fvm->data = ls;

2294:   ls->maxFaces = -1;
2295:   ls->workSize = -1;
2296:   ls->B        = NULL;
2297:   ls->Binv     = NULL;
2298:   ls->tau      = NULL;
2299:   ls->work     = NULL;

2301:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2302:   PetscFVInitialize_LeastSquares(fvm);
2303:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2304:   return(0);
2305: }

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

2312:   Not collective

2314:   Input parameters:
2315: + fvm      - The PetscFV object
2316: - maxFaces - The maximum number of cell faces

2318:   Level: intermediate

2320: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2321: @*/
2322: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2323: {

2328:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2329:   return(0);
2330: }