Actual source code: dtfv.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: #include <petsc-private/petscfvimpl.h> /*I "petscfv.h" I*/
  2: #include <petscdmplex.h>

  4: PetscClassId PETSCLIMITER_CLASSID = 0;

  6: PetscFunctionList PetscLimiterList              = NULL;
  7: PetscBool         PetscLimiterRegisterAllCalled = PETSC_FALSE;

  9: PetscBool Limitercite = PETSC_FALSE;
 10: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
 11:                                "  title   = {Analysis of slope limiters on irregular grids},\n"
 12:                                "  journal = {AIAA paper},\n"
 13:                                "  author  = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
 14:                                "  volume  = {490},\n"
 15:                                "  year    = {2005}\n}\n";

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

 22:   Not Collective

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

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

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

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

 46:   Level: advanced

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

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

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

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

 66:   Collective on PetscLimiter

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

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

 75:   Level: intermediate

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

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

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

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

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

109:   Not Collective

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

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

117:   Level: intermediate

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

129:   if (!PetscLimiterRegisterAllCalled) {PetscLimiterRegisterAll();}
130:   *name = ((PetscObject) lim)->type_name;
131:   return(0);
132: }

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

139:   Collective on PetscLimiter

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

145:   Level: developer

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

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

162: /*
163:   PetscLimiterViewFromOptions - Processes command line options to determine if/how a PetscLimiter is to be viewed.

165:   Collective on PetscLimiter

167:   Input Parameters:
168: + lim    - the PetscLimiter
169: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
170: - optionname - option to activate viewing

172:   Level: intermediate

174: .keywords: PetscLimiter, view, options, database
175: .seealso: VecViewFromOptions(), MatViewFromOptions()
176: */
177: PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter lim, const char prefix[], const char optionname[])
178: {
179:   PetscViewer       viewer;
180:   PetscViewerFormat format;
181:   PetscBool         flg;
182:   PetscErrorCode    ierr;

185:   if (prefix) {PetscOptionsGetViewer(PetscObjectComm((PetscObject) lim), prefix,                      optionname, &viewer, &format, &flg);}
186:   else        {PetscOptionsGetViewer(PetscObjectComm((PetscObject) lim), ((PetscObject) lim)->prefix, optionname, &viewer, &format, &flg);}
187:   if (flg) {
188:     PetscViewerPushFormat(viewer, format);
189:     PetscLimiterView(lim, viewer);
190:     PetscViewerPopFormat(viewer);
191:     PetscViewerDestroy(&viewer);
192:   }
193:   return(0);
194: }

198: /*@
199:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

201:   Collective on PetscLimiter

203:   Input Parameter:
204: . lim - the PetscLimiter object to set options for

206:   Level: developer

208: .seealso: PetscLimiterView()
209: @*/
210: PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
211: {
212:   const char    *defaultType;
213:   char           name[256];
214:   PetscBool      flg;

219:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
220:   else                                 defaultType = ((PetscObject) lim)->type_name;
221:   if (!PetscLimiterRegisterAllCalled) {PetscLimiterRegisterAll();}

223:   PetscObjectOptionsBegin((PetscObject) lim);
224:   PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);
225:   if (flg) {
226:     PetscLimiterSetType(lim, name);
227:   } else if (!((PetscObject) lim)->type_name) {
228:     PetscLimiterSetType(lim, defaultType);
229:   }
230:   if (lim->ops->setfromoptions) {(*lim->ops->setfromoptions)(lim);}
231:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
232:   PetscObjectProcessOptionsHandlers((PetscObject) lim);
233:   PetscOptionsEnd();
234:   PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");
235:   return(0);
236: }

240: /*@C
241:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

243:   Collective on PetscLimiter

245:   Input Parameter:
246: . lim - the PetscLimiter object to setup

248:   Level: developer

250: .seealso: PetscLimiterView(), PetscLimiterDestroy()
251: @*/
252: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
253: {

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

264: /*@
265:   PetscLimiterDestroy - Destroys a PetscLimiter object

267:   Collective on PetscLimiter

269:   Input Parameter:
270: . lim - the PetscLimiter object to destroy

272:   Level: developer

274: .seealso: PetscLimiterView()
275: @*/
276: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
277: {

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

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

287:   if ((*lim)->ops->destroy) {(*(*lim)->ops->destroy)(*lim);}
288:   PetscHeaderDestroy(lim);
289:   return(0);
290: }

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

297:   Collective on MPI_Comm

299:   Input Parameter:
300: . comm - The communicator for the PetscLimiter object

302:   Output Parameter:
303: . lim - The PetscLimiter object

305:   Level: beginner

307: .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
308: @*/
309: PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
310: {
311:   PetscLimiter   l;

316:   PetscCitationsRegister(LimiterCitation,&Limitercite);
317:   *lim = NULL;
318:   PetscFVInitializePackage();

320:   PetscHeaderCreate(l, _p_PetscLimiter, struct _PetscLimiterOps, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);
321:   PetscMemzero(l->ops, sizeof(struct _PetscLimiterOps));

323:   *lim = l;
324:   return(0);
325: }

329: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
330:  *
331:  * The classical flux-limited formulation is psi(r) where
332:  *
333:  * r = (u[0] - u[-1]) / (u[1] - u[0])
334:  *
335:  * The second order TVD region is bounded by
336:  *
337:  * psi_minmod(r) = min(r,1)      and        psi_superbee(r) = min(2, 2r, max(1,r))
338:  *
339:  * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
340:  * phi(r)(r+1)/2 in which the reconstructed interface values are
341:  *
342:  * u(v) = u[0] + phi(r) (grad u)[0] v
343:  *
344:  * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
345:  *
346:  * 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))
347:  *
348:  * For a nicer symmetric formulation, rewrite in terms of
349:  *
350:  * f = (u[0] - u[-1]) / (u[1] - u[-1])
351:  *
352:  * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
353:  *
354:  * phi(r) = phi(1/r)
355:  *
356:  * becomes
357:  *
358:  * w(f) = w(1-f).
359:  *
360:  * The limiters below implement this final form w(f). The reference methods are
361:  *
362:  * w_minmod(f) = 2 min(f,(1-f))             w_superbee(r) = 4 min((1-f), f)
363:  * */
364: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
365: {

371:   (*lim->ops->limit)(lim, flim, phi);
372:   return(0);
373: }

377: PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter fvm)
378: {
379:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) fvm->data;
380:   PetscErrorCode    ierr;

383:   PetscFree(l);
384:   return(0);
385: }

389: PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter fv, PetscViewer viewer)
390: {
391:   PetscViewerFormat format;
392:   PetscErrorCode    ierr;

395:   PetscViewerGetFormat(viewer, &format);
396:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
397:   return(0);
398: }

402: PetscErrorCode PetscLimiterView_Sin(PetscLimiter fv, PetscViewer viewer)
403: {
404:   PetscBool      iascii;

410:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
411:   if (iascii) {PetscLimiterView_Sin_Ascii(fv, viewer);}
412:   return(0);
413: }

417: PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
418: {
420:   *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
421:   return(0);
422: }

426: PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter fvm)
427: {
429:   fvm->ops->view    = PetscLimiterView_Sin;
430:   fvm->ops->destroy = PetscLimiterDestroy_Sin;
431:   fvm->ops->limit   = PetscLimiterLimit_Sin;
432:   return(0);
433: }

435: /*MC
436:   PETSCLIMITERSIN = "sin" - A PetscLimiter object

438:   Level: intermediate

440: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
441: M*/

445: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
446: {
447:   PetscLimiter_Sin *l;
448:   PetscErrorCode    ierr;

452:   PetscNewLog(lim, &l);
453:   lim->data = l;

455:   PetscLimiterInitialize_Sin(lim);
456:   return(0);
457: }

461: PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter fvm)
462: {
463:   PetscLimiter_Zero *l = (PetscLimiter_Zero *) fvm->data;
464:   PetscErrorCode    ierr;

467:   PetscFree(l);
468:   return(0);
469: }

473: PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter fv, PetscViewer viewer)
474: {
475:   PetscViewerFormat format;
476:   PetscErrorCode    ierr;

479:   PetscViewerGetFormat(viewer, &format);
480:   PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
481:   return(0);
482: }

486: PetscErrorCode PetscLimiterView_Zero(PetscLimiter fv, PetscViewer viewer)
487: {
488:   PetscBool      iascii;

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

501: PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
502: {
504:   *phi = 0.0;
505:   return(0);
506: }

510: PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter fvm)
511: {
513:   fvm->ops->view    = PetscLimiterView_Zero;
514:   fvm->ops->destroy = PetscLimiterDestroy_Zero;
515:   fvm->ops->limit   = PetscLimiterLimit_Zero;
516:   return(0);
517: }

519: /*MC
520:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

522:   Level: intermediate

524: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
525: M*/

529: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
530: {
531:   PetscLimiter_Zero *l;
532:   PetscErrorCode     ierr;

536:   PetscNewLog(lim, &l);
537:   lim->data = l;

539:   PetscLimiterInitialize_Zero(lim);
540:   return(0);
541: }

545: PetscErrorCode PetscLimiterDestroy_None(PetscLimiter fvm)
546: {
547:   PetscLimiter_None *l = (PetscLimiter_None *) fvm->data;
548:   PetscErrorCode    ierr;

551:   PetscFree(l);
552:   return(0);
553: }

557: PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter fv, PetscViewer viewer)
558: {
559:   PetscViewerFormat format;
560:   PetscErrorCode    ierr;

563:   PetscViewerGetFormat(viewer, &format);
564:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
565:   return(0);
566: }

570: PetscErrorCode PetscLimiterView_None(PetscLimiter fv, PetscViewer viewer)
571: {
572:   PetscBool      iascii;

578:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
579:   if (iascii) {PetscLimiterView_None_Ascii(fv, viewer);}
580:   return(0);
581: }

585: PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
586: {
588:   *phi = 1.0;
589:   return(0);
590: }

594: PetscErrorCode PetscLimiterInitialize_None(PetscLimiter fvm)
595: {
597:   fvm->ops->view    = PetscLimiterView_None;
598:   fvm->ops->destroy = PetscLimiterDestroy_None;
599:   fvm->ops->limit   = PetscLimiterLimit_None;
600:   return(0);
601: }

603: /*MC
604:   PETSCLIMITERNONE = "none" - A PetscLimiter object

606:   Level: intermediate

608: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
609: M*/

613: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
614: {
615:   PetscLimiter_None *l;
616:   PetscErrorCode    ierr;

620:   PetscNewLog(lim, &l);
621:   lim->data = l;

623:   PetscLimiterInitialize_None(lim);
624:   return(0);
625: }

629: PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter fvm)
630: {
631:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) fvm->data;
632:   PetscErrorCode    ierr;

635:   PetscFree(l);
636:   return(0);
637: }

641: PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter fv, PetscViewer viewer)
642: {
643:   PetscViewerFormat format;
644:   PetscErrorCode    ierr;

647:   PetscViewerGetFormat(viewer, &format);
648:   PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");
649:   return(0);
650: }

654: PetscErrorCode PetscLimiterView_Minmod(PetscLimiter fv, PetscViewer viewer)
655: {
656:   PetscBool      iascii;

662:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
663:   if (iascii) {PetscLimiterView_Minmod_Ascii(fv, viewer);}
664:   return(0);
665: }

669: PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
670: {
672:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
673:   return(0);
674: }

678: PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter fvm)
679: {
681:   fvm->ops->view    = PetscLimiterView_Minmod;
682:   fvm->ops->destroy = PetscLimiterDestroy_Minmod;
683:   fvm->ops->limit   = PetscLimiterLimit_Minmod;
684:   return(0);
685: }

687: /*MC
688:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

690:   Level: intermediate

692: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
693: M*/

697: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
698: {
699:   PetscLimiter_Minmod *l;
700:   PetscErrorCode    ierr;

704:   PetscNewLog(lim, &l);
705:   lim->data = l;

707:   PetscLimiterInitialize_Minmod(lim);
708:   return(0);
709: }

713: PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter fvm)
714: {
715:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) fvm->data;
716:   PetscErrorCode    ierr;

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

725: PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter fv, PetscViewer viewer)
726: {
727:   PetscViewerFormat format;
728:   PetscErrorCode    ierr;

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

738: PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter fv, PetscViewer viewer)
739: {
740:   PetscBool      iascii;

746:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
747:   if (iascii) {PetscLimiterView_VanLeer_Ascii(fv, viewer);}
748:   return(0);
749: }

753: PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
754: {
756:   *phi = PetscMax(0, 4*f*(1-f));
757:   return(0);
758: }

762: PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter fvm)
763: {
765:   fvm->ops->view    = PetscLimiterView_VanLeer;
766:   fvm->ops->destroy = PetscLimiterDestroy_VanLeer;
767:   fvm->ops->limit   = PetscLimiterLimit_VanLeer;
768:   return(0);
769: }

771: /*MC
772:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

774:   Level: intermediate

776: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
777: M*/

781: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
782: {
783:   PetscLimiter_VanLeer *l;
784:   PetscErrorCode    ierr;

788:   PetscNewLog(lim, &l);
789:   lim->data = l;

791:   PetscLimiterInitialize_VanLeer(lim);
792:   return(0);
793: }

797: PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter fvm)
798: {
799:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) fvm->data;
800:   PetscErrorCode    ierr;

803:   PetscFree(l);
804:   return(0);
805: }

809: PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter fv, PetscViewer viewer)
810: {
811:   PetscViewerFormat format;
812:   PetscErrorCode    ierr;

815:   PetscViewerGetFormat(viewer, &format);
816:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
817:   return(0);
818: }

822: PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter fv, PetscViewer viewer)
823: {
824:   PetscBool      iascii;

830:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
831:   if (iascii) {PetscLimiterView_VanAlbada_Ascii(fv, viewer);}
832:   return(0);
833: }

837: PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
838: {
840:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
841:   return(0);
842: }

846: PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter fvm)
847: {
849:   fvm->ops->view    = PetscLimiterView_VanAlbada;
850:   fvm->ops->destroy = PetscLimiterDestroy_VanAlbada;
851:   fvm->ops->limit   = PetscLimiterLimit_VanAlbada;
852:   return(0);
853: }

855: /*MC
856:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

858:   Level: intermediate

860: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
861: M*/

865: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
866: {
867:   PetscLimiter_VanAlbada *l;
868:   PetscErrorCode    ierr;

872:   PetscNewLog(lim, &l);
873:   lim->data = l;

875:   PetscLimiterInitialize_VanAlbada(lim);
876:   return(0);
877: }

881: PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter fvm)
882: {
883:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) fvm->data;
884:   PetscErrorCode    ierr;

887:   PetscFree(l);
888:   return(0);
889: }

893: PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter fv, PetscViewer viewer)
894: {
895:   PetscViewerFormat format;
896:   PetscErrorCode    ierr;

899:   PetscViewerGetFormat(viewer, &format);
900:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
901:   return(0);
902: }

906: PetscErrorCode PetscLimiterView_Superbee(PetscLimiter fv, PetscViewer viewer)
907: {
908:   PetscBool      iascii;

914:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
915:   if (iascii) {PetscLimiterView_Superbee_Ascii(fv, viewer);}
916:   return(0);
917: }

921: PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
922: {
924:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
925:   return(0);
926: }

930: PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter fvm)
931: {
933:   fvm->ops->view    = PetscLimiterView_Superbee;
934:   fvm->ops->destroy = PetscLimiterDestroy_Superbee;
935:   fvm->ops->limit   = PetscLimiterLimit_Superbee;
936:   return(0);
937: }

939: /*MC
940:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

942:   Level: intermediate

944: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
945: M*/

949: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
950: {
951:   PetscLimiter_Superbee *l;
952:   PetscErrorCode    ierr;

956:   PetscNewLog(lim, &l);
957:   lim->data = l;

959:   PetscLimiterInitialize_Superbee(lim);
960:   return(0);
961: }

965: PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter fvm)
966: {
967:   PetscLimiter_MC *l = (PetscLimiter_MC *) fvm->data;
968:   PetscErrorCode    ierr;

971:   PetscFree(l);
972:   return(0);
973: }

977: PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter fv, PetscViewer viewer)
978: {
979:   PetscViewerFormat format;
980:   PetscErrorCode    ierr;

983:   PetscViewerGetFormat(viewer, &format);
984:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
985:   return(0);
986: }

990: PetscErrorCode PetscLimiterView_MC(PetscLimiter fv, PetscViewer viewer)
991: {
992:   PetscBool      iascii;

998:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
999:   if (iascii) {PetscLimiterView_MC_Ascii(fv, viewer);}
1000:   return(0);
1001: }

1005: /* aka Barth-Jespersen */
1006: PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
1007: {
1009:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
1010:   return(0);
1011: }

1015: PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter fvm)
1016: {
1018:   fvm->ops->view    = PetscLimiterView_MC;
1019:   fvm->ops->destroy = PetscLimiterDestroy_MC;
1020:   fvm->ops->limit   = PetscLimiterLimit_MC;
1021:   return(0);
1022: }

1024: /*MC
1025:   PETSCLIMITERMC = "mc" - A PetscLimiter object

1027:   Level: intermediate

1029: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
1030: M*/

1034: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
1035: {
1036:   PetscLimiter_MC *l;
1037:   PetscErrorCode    ierr;

1041:   PetscNewLog(lim, &l);
1042:   lim->data = l;

1044:   PetscLimiterInitialize_MC(lim);
1045:   return(0);
1046: }

1048: PetscClassId PETSCFV_CLASSID = 0;

1050: PetscFunctionList PetscFVList              = NULL;
1051: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

1055: /*@C
1056:   PetscFVRegister - Adds a new PetscFV implementation

1058:   Not Collective

1060:   Input Parameters:
1061: + name        - The name of a new user-defined creation routine
1062: - create_func - The creation routine itself

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

1067:   Sample usage:
1068: .vb
1069:     PetscFVRegister("my_fv", MyPetscFVCreate);
1070: .ve

1072:   Then, your PetscFV type can be chosen with the procedural interface via
1073: .vb
1074:     PetscFVCreate(MPI_Comm, PetscFV *);
1075:     PetscFVSetType(PetscFV, "my_fv");
1076: .ve
1077:    or at runtime via the option
1078: .vb
1079:     -petscfv_type my_fv
1080: .ve

1082:   Level: advanced

1084: .keywords: PetscFV, register
1085: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()

1087: @*/
1088: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
1089: {

1093:   PetscFunctionListAdd(&PetscFVList, sname, function);
1094:   return(0);
1095: }

1099: /*@C
1100:   PetscFVSetType - Builds a particular PetscFV

1102:   Collective on PetscFV

1104:   Input Parameters:
1105: + fvm  - The PetscFV object
1106: - name - The kind of FVM space

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

1111:   Level: intermediate

1113: .keywords: PetscFV, set, type
1114: .seealso: PetscFVGetType(), PetscFVCreate()
1115: @*/
1116: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
1117: {
1118:   PetscErrorCode (*r)(PetscFV);
1119:   PetscBool      match;

1124:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
1125:   if (match) return(0);

1127:   if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}
1128:   PetscFunctionListFind(PetscFVList, name, &r);
1129:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);

1131:   if (fvm->ops->destroy) {
1132:     (*fvm->ops->destroy)(fvm);
1133:     fvm->ops->destroy = NULL;
1134:   }
1135:   (*r)(fvm);
1136:   PetscObjectChangeTypeName((PetscObject) fvm, name);
1137:   return(0);
1138: }

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

1145:   Not Collective

1147:   Input Parameter:
1148: . fvm  - The PetscFV

1150:   Output Parameter:
1151: . name - The PetscFV type name

1153:   Level: intermediate

1155: .keywords: PetscFV, get, type, name
1156: .seealso: PetscFVSetType(), PetscFVCreate()
1157: @*/
1158: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1159: {

1165:   if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}
1166:   *name = ((PetscObject) fvm)->type_name;
1167:   return(0);
1168: }

1172: /*@C
1173:   PetscFVView - Views a PetscFV

1175:   Collective on PetscFV

1177:   Input Parameter:
1178: + fvm - the PetscFV object to view
1179: - v   - the viewer

1181:   Level: developer

1183: .seealso: PetscFVDestroy()
1184: @*/
1185: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1186: {

1191:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);}
1192:   if (fvm->ops->view) {(*fvm->ops->view)(fvm, v);}
1193:   return(0);
1194: }

1198: /*
1199:   PetscFVViewFromOptions - Processes command line options to determine if/how a PetscFV is to be viewed.

1201:   Collective on PetscFV

1203:   Input Parameters:
1204: + fvm    - the PetscFV
1205: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
1206: - optionname - option to activate viewing

1208:   Level: intermediate

1210: .keywords: PetscFV, view, options, database
1211: .seealso: VecViewFromOptions(), MatViewFromOptions()
1212: */
1213: PetscErrorCode PetscFVViewFromOptions(PetscFV fvm, const char prefix[], const char optionname[])
1214: {
1215:   PetscViewer       viewer;
1216:   PetscViewerFormat format;
1217:   PetscBool         flg;
1218:   PetscErrorCode    ierr;

1221:   if (prefix) {PetscOptionsGetViewer(PetscObjectComm((PetscObject) fvm), prefix,                      optionname, &viewer, &format, &flg);}
1222:   else        {PetscOptionsGetViewer(PetscObjectComm((PetscObject) fvm), ((PetscObject) fvm)->prefix, optionname, &viewer, &format, &flg);}
1223:   if (flg) {
1224:     PetscViewerPushFormat(viewer, format);
1225:     PetscFVView(fvm, viewer);
1226:     PetscViewerPopFormat(viewer);
1227:     PetscViewerDestroy(&viewer);
1228:   }
1229:   return(0);
1230: }

1234: /*@
1235:   PetscFVSetFromOptions - sets parameters in a PetscFV from the options database

1237:   Collective on PetscFV

1239:   Input Parameter:
1240: . fvm - the PetscFV object to set options for

1242:   Level: developer

1244: .seealso: PetscFVView()
1245: @*/
1246: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1247: {
1248:   const char    *defaultType;
1249:   char           name[256];
1250:   PetscBool      flg;

1255:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1256:   else                                 defaultType = ((PetscObject) fvm)->type_name;
1257:   if (!PetscFVRegisterAllCalled) {PetscFVRegisterAll();}

1259:   PetscObjectOptionsBegin((PetscObject) fvm);
1260:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1261:   if (flg) {
1262:     PetscFVSetType(fvm, name);
1263:   } else if (!((PetscObject) fvm)->type_name) {
1264:     PetscFVSetType(fvm, defaultType);
1265:   }
1266:   if (fvm->ops->setfromoptions) {(*fvm->ops->setfromoptions)(fvm);}
1267:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1268:   PetscObjectProcessOptionsHandlers((PetscObject) fvm);
1269:   PetscLimiterSetFromOptions(fvm->limiter);
1270:   PetscOptionsEnd();
1271:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1272:   return(0);
1273: }

1277: /*@
1278:   PetscFVSetUp - Construct data structures for the PetscFV

1280:   Collective on PetscFV

1282:   Input Parameter:
1283: . fvm - the PetscFV object to setup

1285:   Level: developer

1287: .seealso: PetscFVView(), PetscFVDestroy()
1288: @*/
1289: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1290: {

1295:   PetscLimiterSetUp(fvm->limiter);
1296:   if (fvm->ops->setup) {(*fvm->ops->setup)(fvm);}
1297:   return(0);
1298: }

1302: /*@
1303:   PetscFVDestroy - Destroys a PetscFV object

1305:   Collective on PetscFV

1307:   Input Parameter:
1308: . fvm - the PetscFV object to destroy

1310:   Level: developer

1312: .seealso: PetscFVView()
1313: @*/
1314: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1315: {

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

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

1325:   PetscLimiterDestroy(&(*fvm)->limiter);
1326:   PetscFree((*fvm)->fluxWork);

1328:   if ((*fvm)->ops->destroy) {(*(*fvm)->ops->destroy)(*fvm);}
1329:   PetscHeaderDestroy(fvm);
1330:   return(0);
1331: }

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

1338:   Collective on MPI_Comm

1340:   Input Parameter:
1341: . comm - The communicator for the PetscFV object

1343:   Output Parameter:
1344: . fvm - The PetscFV object

1346:   Level: beginner

1348: .seealso: PetscFVSetType(), PETSCFVUPWIND
1349: @*/
1350: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1351: {
1352:   PetscFV        f;

1357:   *fvm = NULL;
1358:   PetscFVInitializePackage();

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

1363:   PetscLimiterCreate(comm, &f->limiter);
1364:   f->numComponents    = 1;
1365:   f->dim              = 0;
1366:   f->computeGradients = PETSC_FALSE;
1367:   f->fluxWork         = NULL;

1369:   *fvm = f;
1370:   return(0);
1371: }

1375: /*@
1376:   PetscFVSetLimiter - Set the limiter object

1378:   Logically collective on PetscFV

1380:   Input Parameters:
1381: + fvm - the PetscFV object to destroy
1382: - lim - The PetscLimiter

1384:   Level: developer

1386: .seealso: PetscFVGetLimiter()
1387: @*/
1388: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1389: {

1395:   PetscLimiterDestroy(&fvm->limiter);
1396:   PetscObjectReference((PetscObject) lim);
1397:   fvm->limiter = lim;
1398:   return(0);
1399: }

1403: /*@
1404:   PetscFVGetLimiter - Get the limiter object

1406:   Not collective

1408:   Input Parameter:
1409: . fvm - the PetscFV object to destroy

1411:   Output Parameter:
1412: . lim - The PetscLimiter

1414:   Level: developer

1416: .seealso: PetscFVSetLimiter()
1417: @*/
1418: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1419: {
1423:   *lim = fvm->limiter;
1424:   return(0);
1425: }

1429: /*@
1430:   PetscFVSetNumComponents - Set the number of field components

1432:   Logically collective on PetscFV

1434:   Input Parameters:
1435: + fvm - the PetscFV object to destroy
1436: - comp - The number of components

1438:   Level: developer

1440: .seealso: PetscFVGetNumComponents()
1441: @*/
1442: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1443: {

1448:   fvm->numComponents = comp;
1449:   PetscFree(fvm->fluxWork);
1450:   PetscMalloc1(comp, &fvm->fluxWork);
1451:   return(0);
1452: }

1456: /*@
1457:   PetscFVGetNumComponents - Get the number of field components

1459:   Not collective

1461:   Input Parameter:
1462: . fvm - the PetscFV object to destroy

1464:   Output Parameter:
1465: , comp - The number of components

1467:   Level: developer

1469: .seealso: PetscFVSetNumComponents()
1470: @*/
1471: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1472: {
1476:   *comp = fvm->numComponents;
1477:   return(0);
1478: }

1482: /*@
1483:   PetscFVSetSpatialDimension - Set the spatial dimension

1485:   Logically collective on PetscFV

1487:   Input Parameters:
1488: + fvm - the PetscFV object to destroy
1489: - dim - The spatial dimension

1491:   Level: developer

1493: .seealso: PetscFVGetSpatialDimension()
1494: @*/
1495: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1496: {
1499:   fvm->dim = dim;
1500:   return(0);
1501: }

1505: /*@
1506:   PetscFVGetSpatialDimension - Get the spatial dimension

1508:   Logically collective on PetscFV

1510:   Input Parameter:
1511: . fvm - the PetscFV object to destroy

1513:   Output Parameter:
1514: . dim - The spatial dimension

1516:   Level: developer

1518: .seealso: PetscFVSetSpatialDimension()
1519: @*/
1520: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1521: {
1525:   *dim = fvm->dim;
1526:   return(0);
1527: }

1531: /*@
1532:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1534:   Logically collective on PetscFV

1536:   Input Parameters:
1537: + fvm - the PetscFV object to destroy
1538: - computeGradients - Flag to compute cell gradients

1540:   Level: developer

1542: .seealso: PetscFVGetComputeGradients()
1543: @*/
1544: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1545: {
1548:   fvm->computeGradients = computeGradients;
1549:   return(0);
1550: }

1554: /*@
1555:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1557:   Not collective

1559:   Input Parameter:
1560: . fvm - the PetscFV object to destroy

1562:   Output Parameter:
1563: . computeGradients - Flag to compute cell gradients

1565:   Level: developer

1567: .seealso: PetscFVSetComputeGradients()
1568: @*/
1569: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1570: {
1574:   *computeGradients = fvm->computeGradients;
1575:   return(0);
1576: }

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

1583:   Input Parameters:
1584: + fvm      - The PetscFV object
1585: . numFaces - The number of cell faces which are not constrained
1586: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1588:   Level: developer

1590: .seealso: PetscFVCreate()
1591: @*/
1592: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1593: {

1598:   if (fvm->ops->computegradient) {(*fvm->ops->computegradient)(fvm, numFaces, dx, grad);}
1599:   return(0);
1600: }

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

1607:   Not collective

1609:   Input Parameters:
1610: + fvm          - The PetscFV object for the field being integrated
1611: . Nface        - The number of faces in the chunk
1612: . Nf           - The number of physical fields
1613: . fv           - The PetscFV objects for each field
1614: . field        - The field being integrated
1615: . fgeom        - The face geometry for each face in the chunk
1616: . cgeom        - The cell geometry for each pair of cells in the chunk
1617: . uL           - The state from the cell on the left
1618: . uR           - The state from the cell on the right
1619: . riemann      - Riemann solver
1620: - ctx          - User context passed to Riemann solve

1622:   Output Parameter
1623: + fluxL        - the left fluxes for each face
1624: - fluxR        - the right fluxes for each face
1625: */
1626: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom, PetscScalar uL[], PetscScalar uR[],
1627:                                            void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
1628:                                            PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
1629: {

1634:   if (fvm->ops->integraterhsfunction) {(*fvm->ops->integraterhsfunction)(fvm, Nfaces, Nf, fv, field, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, ctx);}
1635:   return(0);
1636: }

1640: PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1641: {
1642:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1646:   PetscFree(b);
1647:   return(0);
1648: }

1652: PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1653: {
1654:   PetscInt          Nc;
1655:   PetscViewerFormat format;
1656:   PetscErrorCode    ierr;

1659:   PetscFVGetNumComponents(fv, &Nc);
1660:   PetscViewerGetFormat(viewer, &format);
1661:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1662:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1663:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1664:   } else {
1665:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1666:   }
1667:   return(0);
1668: }

1672: PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1673: {
1674:   PetscBool      iascii;

1680:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1681:   if (iascii) {PetscFVView_Upwind_Ascii(fv, viewer);}
1682:   return(0);
1683: }

1687: PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1688: {
1690:   return(0);
1691: }

1695: /*
1696:   fgeom[f]=>v0 is the centroid
1697:   cgeom->vol[f*2+0] contains the left  geom
1698:   cgeom->vol[f*2+1] contains the right geom
1699: */
1700: PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom,
1701:                                                   PetscScalar xL[], PetscScalar xR[],
1702:                                                   void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
1703:                                                   PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
1704: {
1705:   PetscScalar   *flux = fvm->fluxWork;
1706:   PetscInt       dim, pdim, f, d;

1710:   PetscFVGetSpatialDimension(fvm, &dim);
1711:   PetscFVGetNumComponents(fvm, &pdim);
1712:   for (f = 0; f < Nfaces; ++f) {
1713:     (*riemann)(&fgeom.v0[f*dim], &fgeom.n[f*dim], &xL[f*pdim], &xR[f*pdim], flux, ctx);
1714:     for (d = 0; d < pdim; ++d) {
1715:       fluxL[f*pdim+d] = flux[d] / cgeom.vol[f*2+0];
1716:       fluxR[f*pdim+d] = flux[d] / cgeom.vol[f*2+1];
1717:     }
1718:   }
1719:   return(0);
1720: }

1724: PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1725: {
1727:   fvm->ops->setfromoptions          = NULL;
1728:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1729:   fvm->ops->view                    = PetscFVView_Upwind;
1730:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
1731:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
1732:   return(0);
1733: }

1735: /*MC
1736:   PETSCFVUPWIND = "upwind" - A PetscFV object

1738:   Level: intermediate

1740: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1741: M*/

1745: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1746: {
1747:   PetscFV_Upwind *b;
1748:   PetscErrorCode  ierr;

1752:   PetscNewLog(fvm,&b);
1753:   fvm->data = b;

1755:   PetscFVInitialize_Upwind(fvm);
1756:   return(0);
1757: }

1759: #include <petscblaslapack.h>

1763: PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1764: {
1765:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1766:   PetscErrorCode        ierr;

1769:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1770:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1771:   PetscFree(ls);
1772:   return(0);
1773: }

1777: PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1778: {
1779:   PetscInt          Nc;
1780:   PetscViewerFormat format;
1781:   PetscErrorCode    ierr;

1784:   PetscFVGetNumComponents(fv, &Nc);
1785:   PetscViewerGetFormat(viewer, &format);
1786:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1787:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1788:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1789:   } else {
1790:     PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1791:   }
1792:   return(0);
1793: }

1797: PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1798: {
1799:   PetscBool      iascii;

1805:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1806:   if (iascii) {PetscFVView_LeastSquares_Ascii(fv, viewer);}
1807:   return(0);
1808: }

1812: PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1813: {
1815:   return(0);
1816: }

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

1829:   if (debug) {
1830:     PetscMalloc1(m*n,&Aback);
1831:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1832:   }

1834:   PetscBLASIntCast(m,&M);
1835:   PetscBLASIntCast(n,&N);
1836:   PetscBLASIntCast(mstride,&lda);
1837:   PetscBLASIntCast(worksize,&ldwork);
1838:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1839:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1840:   PetscFPTrapPop();
1841:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1842:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1844:   /* Extract an explicit representation of Q */
1845:   Q    = Ainv;
1846:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1847:   K    = N;                     /* full rank */
1848:   LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
1849:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1851:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1852:   Alpha = 1.0;
1853:   ldb   = lda;
1854:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1855:   /* Ainv is Q, overwritten with inverse */

1857:   if (debug) {                      /* Check that pseudo-inverse worked */
1858:     PetscScalar  Beta = 0.0;
1859:     PetscBLASInt ldc;
1860:     K   = N;
1861:     ldc = N;
1862:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1863:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1864:     PetscFree(Aback);
1865:   }
1866:   return(0);
1867: }

1871: /* Overwrites A. Can handle degenerate problems and m<n. */
1872: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1873: {
1874:   PetscBool      debug = PETSC_FALSE;
1875:   PetscScalar   *Brhs, *Aback;
1876:   PetscScalar   *tmpwork;
1877:   PetscReal      rcond;
1878:   PetscInt       i, j, maxmn;
1879:   PetscBLASInt   M, N, nrhs, lda, ldb, irank, ldwork, info;

1883:   if (debug) {
1884:     PetscMalloc1(m*n,&Aback);
1885:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1886:   }

1888:   /* initialize to identity */
1889:   tmpwork = Ainv;
1890:   Brhs = work;
1891:   maxmn = PetscMax(m,n);
1892:   for (j=0; j<maxmn; j++) {
1893:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1894:   }

1896:   PetscBLASIntCast(m,&M);
1897:   PetscBLASIntCast(n,&N);
1898:   nrhs  = M;
1899:   PetscBLASIntCast(mstride,&lda);
1900:   PetscBLASIntCast(maxmn,&ldb);
1901:   PetscBLASIntCast(worksize,&ldwork);
1902:   rcond = -1;
1903:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1904: #if defined(PETSC_USE_COMPLEX)
1905:   if (tmpwork && rcond) rcond = 0.0; /* Get rid of compiler warning */
1906:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "I don't think this makes sense for complex numbers");
1907: #else
1908:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
1909: #endif
1910:   PetscFPTrapPop();
1911:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
1912:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1913:   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");

1915:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
1916:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
1917:   for (i=0; i<n; i++) {
1918:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
1919:   }
1920:   return(0);
1921: }

1923: #if 0
1926: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1927: {
1928:   PetscReal       grad[2] = {0, 0};
1929:   const PetscInt *faces;
1930:   PetscInt        numFaces, f;
1931:   PetscErrorCode  ierr;

1934:   DMPlexGetConeSize(dm, cell, &numFaces);
1935:   DMPlexGetCone(dm, cell, &faces);
1936:   for (f = 0; f < numFaces; ++f) {
1937:     const PetscInt *fcells;
1938:     const CellGeom *cg1;
1939:     const FaceGeom *fg;

1941:     DMPlexGetSupport(dm, faces[f], &fcells);
1942:     DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);
1943:     for (i = 0; i < 2; ++i) {
1944:       PetscScalar du;

1946:       if (fcells[i] == c) continue;
1947:       DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);
1948:       du   = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1949:       grad[0] += fg->grad[!i][0] * du;
1950:       grad[1] += fg->grad[!i][1] * du;
1951:     }
1952:   }
1953:   PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);
1954:   return(0);
1955: }
1956: #endif

1960: /*
1961:   PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell

1963:   Input Parameters:
1964: + fvm      - The PetscFV object
1965: . numFaces - The number of cell faces which are not constrained
1966: . dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1968:   Level: developer

1970: .seealso: PetscFVCreate()
1971: */
1972: PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
1973: {
1974:   PetscFV_LeastSquares *ls       = (PetscFV_LeastSquares *) fvm->data;
1975:   const PetscBool       useSVD   = PETSC_TRUE;
1976:   const PetscInt        maxFaces = ls->maxFaces;
1977:   PetscInt              dim, f, d;
1978:   PetscErrorCode        ierr;

1981:   if (numFaces > maxFaces) {
1982:     if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
1983:     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
1984:   }
1985:   PetscFVGetSpatialDimension(fvm, &dim);
1986:   for (f = 0; f < numFaces; ++f) {
1987:     for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
1988:   }
1989:   /* Overwrites B with garbage, returns Binv in row-major format */
1990:   if (useSVD) {PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
1991:   else        {PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);}
1992:   for (f = 0; f < numFaces; ++f) {
1993:     for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces+f];
1994:   }
1995:   return(0);
1996: }

2000: /*
2001:   fgeom[f]=>v0 is the centroid
2002:   cgeom->vol[f*2+0] contains the left  geom
2003:   cgeom->vol[f*2+1] contains the right geom
2004: */
2005: PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscInt Nfaces, PetscInt Nf, PetscFV fv[], PetscInt field, PetscCellGeometry fgeom, PetscCellGeometry cgeom,
2006:                                                         PetscScalar xL[], PetscScalar xR[],
2007:                                                         void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx),
2008:                                                         PetscScalar fluxL[], PetscScalar fluxR[], void *ctx)
2009: {
2010:   PetscScalar   *flux = fvm->fluxWork;
2011:   PetscInt       dim, pdim, f, d;

2015:   PetscFVGetSpatialDimension(fvm, &dim);
2016:   PetscFVGetNumComponents(fvm, &pdim);
2017:   for (f = 0; f < Nfaces; ++f) {
2018:     (*riemann)(&fgeom.v0[f*dim], &fgeom.n[f*dim], &xL[f*pdim], &xR[f*pdim], flux, ctx);
2019:     for (d = 0; d < pdim; ++d) {
2020:       fluxL[f*pdim+d] = flux[d] / cgeom.vol[f*2+0];
2021:       fluxR[f*pdim+d] = flux[d] / cgeom.vol[f*2+1];
2022:     }
2023:   }
2024:   return(0);
2025: }

2029: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2030: {
2031:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2032:   PetscInt              dim, m, n, nrhs, minwork;
2033:   PetscErrorCode        ierr;

2037:   PetscFVGetSpatialDimension(fvm, &dim);
2038:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2039:   ls->maxFaces = maxFaces;
2040:   m       = ls->maxFaces;
2041:   n       = dim;
2042:   nrhs    = ls->maxFaces;
2043:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
2044:   ls->workSize = 5*minwork; /* We can afford to be extra generous */
2045:   PetscMalloc4(ls->maxFaces*dim,&ls->B,ls->workSize,&ls->Binv,ls->maxFaces,&ls->tau,ls->workSize,&ls->work);
2046:   return(0);
2047: }

2051: PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2052: {
2054:   fvm->ops->setfromoptions          = NULL;
2055:   fvm->ops->setup                   = PetscFVSetUp_LeastSquares;
2056:   fvm->ops->view                    = PetscFVView_LeastSquares;
2057:   fvm->ops->destroy                 = PetscFVDestroy_LeastSquares;
2058:   fvm->ops->computegradient         = PetscFVComputeGradient_LeastSquares;
2059:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_LeastSquares;
2060:   return(0);
2061: }

2063: /*MC
2064:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2066:   Level: intermediate

2068: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2069: M*/

2073: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2074: {
2075:   PetscFV_LeastSquares *ls;
2076:   PetscErrorCode        ierr;

2080:   PetscNewLog(fvm, &ls);
2081:   fvm->data = ls;

2083:   ls->maxFaces = -1;
2084:   ls->workSize = -1;
2085:   ls->B        = NULL;
2086:   ls->Binv     = NULL;
2087:   ls->tau      = NULL;
2088:   ls->work     = NULL;

2090:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
2091:   PetscFVInitialize_LeastSquares(fvm);
2092:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);
2093:   return(0);
2094: }

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

2101:   Not collective

2103:   Input parameters:
2104: + fvm      - The PetscFV object
2105: - maxFaces - The maximum number of cell faces

2107:   Level: intermediate

2109: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2110: @*/
2111: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2112: {

2117:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2118:   return(0);
2119: }