Actual source code: fv.c

  1: #include <petsc/private/petscfvimpl.h>
  2: #include <petscdmplex.h>
  3: #include <petscdmplextransform.h>
  4: #include <petscds.h>

  6: PetscClassId PETSCLIMITER_CLASSID = 0;

  8: PetscFunctionList PetscLimiterList              = NULL;
  9: PetscBool         PetscLimiterRegisterAllCalled = PETSC_FALSE;

 11: PetscBool Limitercite = PETSC_FALSE;
 12: const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
 13:                                "  title   = {Analysis of slope limiters on irregular grids},\n"
 14:                                "  journal = {AIAA paper},\n"
 15:                                "  author  = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
 16:                                "  volume  = {490},\n"
 17:                                "  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: .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()

 50: @*/
 51: PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
 52: {
 53:   PetscFunctionListAdd(&PetscLimiterList, sname, function);
 54:   return 0;
 55: }

 57: /*@C
 58:   PetscLimiterSetType - Builds a particular PetscLimiter

 60:   Collective on lim

 62:   Input Parameters:
 63: + lim  - The PetscLimiter object
 64: - name - The kind of limiter

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

 69:   Level: intermediate

 71: .seealso: PetscLimiterGetType(), PetscLimiterCreate()
 72: @*/
 73: PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
 74: {
 75:   PetscErrorCode (*r)(PetscLimiter);
 76:   PetscBool      match;

 79:   PetscObjectTypeCompare((PetscObject) lim, name, &match);
 80:   if (match) return 0;

 82:   PetscLimiterRegisterAll();
 83:   PetscFunctionListFind(PetscLimiterList, name, &r);

 86:   if (lim->ops->destroy) {
 87:     (*lim->ops->destroy)(lim);
 88:     lim->ops->destroy = NULL;
 89:   }
 90:   (*r)(lim);
 91:   PetscObjectChangeTypeName((PetscObject) lim, name);
 92:   return 0;
 93: }

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

 98:   Not Collective

100:   Input Parameter:
101: . lim  - The PetscLimiter

103:   Output Parameter:
104: . name - The PetscLimiter type name

106:   Level: intermediate

108: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
109: @*/
110: PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
111: {
114:   PetscLimiterRegisterAll();
115:   *name = ((PetscObject) lim)->type_name;
116:   return 0;
117: }

119: /*@C
120:    PetscLimiterViewFromOptions - View from Options

122:    Collective on PetscLimiter

124:    Input Parameters:
125: +  A - the PetscLimiter object to view
126: .  obj - Optional object
127: -  name - command line option

129:    Level: intermediate
130: .seealso:  PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
131: @*/
132: PetscErrorCode  PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
133: {
135:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
136:   return 0;
137: }

139: /*@C
140:   PetscLimiterView - Views a PetscLimiter

142:   Collective on lim

144:   Input Parameters:
145: + lim - the PetscLimiter object to view
146: - v   - the viewer

148:   Level: beginner

150: .seealso: PetscLimiterDestroy()
151: @*/
152: PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
153: {
155:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);
156:   if (lim->ops->view) (*lim->ops->view)(lim, v);
157:   return 0;
158: }

160: /*@
161:   PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database

163:   Collective on lim

165:   Input Parameter:
166: . lim - the PetscLimiter object to set options for

168:   Level: intermediate

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

180:   if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
181:   else                                 defaultType = ((PetscObject) lim)->type_name;
182:   PetscLimiterRegisterAll();

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

199: /*@C
200:   PetscLimiterSetUp - Construct data structures for the PetscLimiter

202:   Collective on lim

204:   Input Parameter:
205: . lim - the PetscLimiter object to setup

207:   Level: intermediate

209: .seealso: PetscLimiterView(), PetscLimiterDestroy()
210: @*/
211: PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
212: {
214:   if (lim->ops->setup) (*lim->ops->setup)(lim);
215:   return 0;
216: }

218: /*@
219:   PetscLimiterDestroy - Destroys a PetscLimiter object

221:   Collective on lim

223:   Input Parameter:
224: . lim - the PetscLimiter object to destroy

226:   Level: beginner

228: .seealso: PetscLimiterView()
229: @*/
230: PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
231: {
232:   if (!*lim) return 0;

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

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

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

246:   Collective

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

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

254:   Level: beginner

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

263:   PetscCitationsRegister(LimiterCitation,&Limitercite);
264:   *lim = NULL;
265:   PetscFVInitializePackage();

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

269:   *lim = l;
270:   return 0;
271: }

273: /*@
274:   PetscLimiterLimit - Limit the flux

276:   Input Parameters:
277: + lim  - The PetscLimiter
278: - flim - The input field

280:   Output Parameter:
281: . phi  - The limited field

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

317:   Level: beginner

319: .seealso: PetscLimiterSetType(), PetscLimiterCreate()
320: @*/
321: PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
322: {
325:   (*lim->ops->limit)(lim, flim, phi);
326:   return 0;
327: }

329: static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
330: {
331:   PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;

333:   PetscFree(l);
334:   return 0;
335: }

337: static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
338: {
339:   PetscViewerFormat format;

341:   PetscViewerGetFormat(viewer, &format);
342:   PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");
343:   return 0;
344: }

346: static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
347: {
348:   PetscBool      iascii;

352:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
353:   if (iascii) PetscLimiterView_Sin_Ascii(lim, viewer);
354:   return 0;
355: }

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

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

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

374:   Level: intermediate

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

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

384:   PetscNewLog(lim, &l);
385:   lim->data = l;

387:   PetscLimiterInitialize_Sin(lim);
388:   return 0;
389: }

391: static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
392: {
393:   PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;

395:   PetscFree(l);
396:   return 0;
397: }

399: static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
400: {
401:   PetscViewerFormat format;

403:   PetscViewerGetFormat(viewer, &format);
404:   PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");
405:   return 0;
406: }

408: static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
409: {
410:   PetscBool      iascii;

414:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
415:   if (iascii) PetscLimiterView_Zero_Ascii(lim, viewer);
416:   return 0;
417: }

419: static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
420: {
421:   *phi = 0.0;
422:   return 0;
423: }

425: static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
426: {
427:   lim->ops->view    = PetscLimiterView_Zero;
428:   lim->ops->destroy = PetscLimiterDestroy_Zero;
429:   lim->ops->limit   = PetscLimiterLimit_Zero;
430:   return 0;
431: }

433: /*MC
434:   PETSCLIMITERZERO = "zero" - A PetscLimiter object

436:   Level: intermediate

438: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
439: M*/

441: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
442: {
443:   PetscLimiter_Zero *l;

446:   PetscNewLog(lim, &l);
447:   lim->data = l;

449:   PetscLimiterInitialize_Zero(lim);
450:   return 0;
451: }

453: static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
454: {
455:   PetscLimiter_None *l = (PetscLimiter_None *) lim->data;

457:   PetscFree(l);
458:   return 0;
459: }

461: static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
462: {
463:   PetscViewerFormat format;

465:   PetscViewerGetFormat(viewer, &format);
466:   PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");
467:   return 0;
468: }

470: static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
471: {
472:   PetscBool      iascii;

476:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
477:   if (iascii) PetscLimiterView_None_Ascii(lim, viewer);
478:   return 0;
479: }

481: static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
482: {
483:   *phi = 1.0;
484:   return 0;
485: }

487: static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
488: {
489:   lim->ops->view    = PetscLimiterView_None;
490:   lim->ops->destroy = PetscLimiterDestroy_None;
491:   lim->ops->limit   = PetscLimiterLimit_None;
492:   return 0;
493: }

495: /*MC
496:   PETSCLIMITERNONE = "none" - A PetscLimiter object

498:   Level: intermediate

500: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
501: M*/

503: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
504: {
505:   PetscLimiter_None *l;

508:   PetscNewLog(lim, &l);
509:   lim->data = l;

511:   PetscLimiterInitialize_None(lim);
512:   return 0;
513: }

515: static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
516: {
517:   PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;

519:   PetscFree(l);
520:   return 0;
521: }

523: static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
524: {
525:   PetscViewerFormat format;

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

532: static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
533: {
534:   PetscBool      iascii;

538:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
539:   if (iascii) PetscLimiterView_Minmod_Ascii(lim, viewer);
540:   return 0;
541: }

543: static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
544: {
545:   *phi = 2*PetscMax(0, PetscMin(f, 1-f));
546:   return 0;
547: }

549: static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
550: {
551:   lim->ops->view    = PetscLimiterView_Minmod;
552:   lim->ops->destroy = PetscLimiterDestroy_Minmod;
553:   lim->ops->limit   = PetscLimiterLimit_Minmod;
554:   return 0;
555: }

557: /*MC
558:   PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object

560:   Level: intermediate

562: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
563: M*/

565: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
566: {
567:   PetscLimiter_Minmod *l;

570:   PetscNewLog(lim, &l);
571:   lim->data = l;

573:   PetscLimiterInitialize_Minmod(lim);
574:   return 0;
575: }

577: static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
578: {
579:   PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;

581:   PetscFree(l);
582:   return 0;
583: }

585: static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
586: {
587:   PetscViewerFormat format;

589:   PetscViewerGetFormat(viewer, &format);
590:   PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");
591:   return 0;
592: }

594: static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
595: {
596:   PetscBool      iascii;

600:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
601:   if (iascii) PetscLimiterView_VanLeer_Ascii(lim, viewer);
602:   return 0;
603: }

605: static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
606: {
607:   *phi = PetscMax(0, 4*f*(1-f));
608:   return 0;
609: }

611: static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
612: {
613:   lim->ops->view    = PetscLimiterView_VanLeer;
614:   lim->ops->destroy = PetscLimiterDestroy_VanLeer;
615:   lim->ops->limit   = PetscLimiterLimit_VanLeer;
616:   return 0;
617: }

619: /*MC
620:   PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object

622:   Level: intermediate

624: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
625: M*/

627: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
628: {
629:   PetscLimiter_VanLeer *l;

632:   PetscNewLog(lim, &l);
633:   lim->data = l;

635:   PetscLimiterInitialize_VanLeer(lim);
636:   return 0;
637: }

639: static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
640: {
641:   PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;

643:   PetscFree(l);
644:   return 0;
645: }

647: static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
648: {
649:   PetscViewerFormat format;

651:   PetscViewerGetFormat(viewer, &format);
652:   PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");
653:   return 0;
654: }

656: static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
657: {
658:   PetscBool      iascii;

662:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
663:   if (iascii) PetscLimiterView_VanAlbada_Ascii(lim, viewer);
664:   return 0;
665: }

667: static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
668: {
669:   *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
670:   return 0;
671: }

673: static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
674: {
675:   lim->ops->view    = PetscLimiterView_VanAlbada;
676:   lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
677:   lim->ops->limit   = PetscLimiterLimit_VanAlbada;
678:   return 0;
679: }

681: /*MC
682:   PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object

684:   Level: intermediate

686: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
687: M*/

689: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
690: {
691:   PetscLimiter_VanAlbada *l;

694:   PetscNewLog(lim, &l);
695:   lim->data = l;

697:   PetscLimiterInitialize_VanAlbada(lim);
698:   return 0;
699: }

701: static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
702: {
703:   PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;

705:   PetscFree(l);
706:   return 0;
707: }

709: static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
710: {
711:   PetscViewerFormat format;

713:   PetscViewerGetFormat(viewer, &format);
714:   PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");
715:   return 0;
716: }

718: static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
719: {
720:   PetscBool      iascii;

724:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
725:   if (iascii) PetscLimiterView_Superbee_Ascii(lim, viewer);
726:   return 0;
727: }

729: static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
730: {
731:   *phi = 4*PetscMax(0, PetscMin(f, 1-f));
732:   return 0;
733: }

735: static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
736: {
737:   lim->ops->view    = PetscLimiterView_Superbee;
738:   lim->ops->destroy = PetscLimiterDestroy_Superbee;
739:   lim->ops->limit   = PetscLimiterLimit_Superbee;
740:   return 0;
741: }

743: /*MC
744:   PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object

746:   Level: intermediate

748: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
749: M*/

751: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
752: {
753:   PetscLimiter_Superbee *l;

756:   PetscNewLog(lim, &l);
757:   lim->data = l;

759:   PetscLimiterInitialize_Superbee(lim);
760:   return 0;
761: }

763: static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
764: {
765:   PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;

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

771: static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
772: {
773:   PetscViewerFormat format;

775:   PetscViewerGetFormat(viewer, &format);
776:   PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");
777:   return 0;
778: }

780: static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
781: {
782:   PetscBool      iascii;

786:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
787:   if (iascii) PetscLimiterView_MC_Ascii(lim, viewer);
788:   return 0;
789: }

791: /* aka Barth-Jespersen */
792: static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
793: {
794:   *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
795:   return 0;
796: }

798: static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
799: {
800:   lim->ops->view    = PetscLimiterView_MC;
801:   lim->ops->destroy = PetscLimiterDestroy_MC;
802:   lim->ops->limit   = PetscLimiterLimit_MC;
803:   return 0;
804: }

806: /*MC
807:   PETSCLIMITERMC = "mc" - A PetscLimiter object

809:   Level: intermediate

811: .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
812: M*/

814: PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
815: {
816:   PetscLimiter_MC *l;

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

822:   PetscLimiterInitialize_MC(lim);
823:   return 0;
824: }

826: PetscClassId PETSCFV_CLASSID = 0;

828: PetscFunctionList PetscFVList              = NULL;
829: PetscBool         PetscFVRegisterAllCalled = PETSC_FALSE;

831: /*@C
832:   PetscFVRegister - Adds a new PetscFV implementation

834:   Not Collective

836:   Input Parameters:
837: + name        - The name of a new user-defined creation routine
838: - create_func - The creation routine itself

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

843:   Sample usage:
844: .vb
845:     PetscFVRegister("my_fv", MyPetscFVCreate);
846: .ve

848:   Then, your PetscFV type can be chosen with the procedural interface via
849: .vb
850:     PetscFVCreate(MPI_Comm, PetscFV *);
851:     PetscFVSetType(PetscFV, "my_fv");
852: .ve
853:    or at runtime via the option
854: .vb
855:     -petscfv_type my_fv
856: .ve

858:   Level: advanced

860: .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()

862: @*/
863: PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
864: {
865:   PetscFunctionListAdd(&PetscFVList, sname, function);
866:   return 0;
867: }

869: /*@C
870:   PetscFVSetType - Builds a particular PetscFV

872:   Collective on fvm

874:   Input Parameters:
875: + fvm  - The PetscFV object
876: - name - The kind of FVM space

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

881:   Level: intermediate

883: .seealso: PetscFVGetType(), PetscFVCreate()
884: @*/
885: PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
886: {
887:   PetscErrorCode (*r)(PetscFV);
888:   PetscBool      match;

891:   PetscObjectTypeCompare((PetscObject) fvm, name, &match);
892:   if (match) return 0;

894:   PetscFVRegisterAll();
895:   PetscFunctionListFind(PetscFVList, name, &r);

898:   if (fvm->ops->destroy) {
899:     (*fvm->ops->destroy)(fvm);
900:     fvm->ops->destroy = NULL;
901:   }
902:   (*r)(fvm);
903:   PetscObjectChangeTypeName((PetscObject) fvm, name);
904:   return 0;
905: }

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

910:   Not Collective

912:   Input Parameter:
913: . fvm  - The PetscFV

915:   Output Parameter:
916: . name - The PetscFV type name

918:   Level: intermediate

920: .seealso: PetscFVSetType(), PetscFVCreate()
921: @*/
922: PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
923: {
926:   PetscFVRegisterAll();
927:   *name = ((PetscObject) fvm)->type_name;
928:   return 0;
929: }

931: /*@C
932:    PetscFVViewFromOptions - View from Options

934:    Collective on PetscFV

936:    Input Parameters:
937: +  A - the PetscFV object
938: .  obj - Optional object
939: -  name - command line option

941:    Level: intermediate
942: .seealso:  PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
943: @*/
944: PetscErrorCode  PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
945: {
947:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
948:   return 0;
949: }

951: /*@C
952:   PetscFVView - Views a PetscFV

954:   Collective on fvm

956:   Input Parameters:
957: + fvm - the PetscFV object to view
958: - v   - the viewer

960:   Level: beginner

962: .seealso: PetscFVDestroy()
963: @*/
964: PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
965: {
967:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);
968:   if (fvm->ops->view) (*fvm->ops->view)(fvm, v);
969:   return 0;
970: }

972: /*@
973:   PetscFVSetFromOptions - sets parameters in a PetscFV from the options database

975:   Collective on fvm

977:   Input Parameter:
978: . fvm - the PetscFV object to set options for

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

983:   Level: intermediate

985: .seealso: PetscFVView()
986: @*/
987: PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
988: {
989:   const char    *defaultType;
990:   char           name[256];
991:   PetscBool      flg;

995:   if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
996:   else                                 defaultType = ((PetscObject) fvm)->type_name;
997:   PetscFVRegisterAll();

999:   PetscObjectOptionsBegin((PetscObject) fvm);
1000:   PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);
1001:   if (flg) {
1002:     PetscFVSetType(fvm, name);
1003:   } else if (!((PetscObject) fvm)->type_name) {
1004:     PetscFVSetType(fvm, defaultType);

1006:   }
1007:   PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);
1008:   if (fvm->ops->setfromoptions) (*fvm->ops->setfromoptions)(fvm);
1009:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1010:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);
1011:   PetscLimiterSetFromOptions(fvm->limiter);
1012:   PetscOptionsEnd();
1013:   PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");
1014:   return 0;
1015: }

1017: /*@
1018:   PetscFVSetUp - Construct data structures for the PetscFV

1020:   Collective on fvm

1022:   Input Parameter:
1023: . fvm - the PetscFV object to setup

1025:   Level: intermediate

1027: .seealso: PetscFVView(), PetscFVDestroy()
1028: @*/
1029: PetscErrorCode PetscFVSetUp(PetscFV fvm)
1030: {
1032:   PetscLimiterSetUp(fvm->limiter);
1033:   if (fvm->ops->setup) (*fvm->ops->setup)(fvm);
1034:   return 0;
1035: }

1037: /*@
1038:   PetscFVDestroy - Destroys a PetscFV object

1040:   Collective on fvm

1042:   Input Parameter:
1043: . fvm - the PetscFV object to destroy

1045:   Level: beginner

1047: .seealso: PetscFVView()
1048: @*/
1049: PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1050: {
1051:   PetscInt       i;

1053:   if (!*fvm) return 0;

1056:   if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; return 0;}
1057:   ((PetscObject) (*fvm))->refct = 0;

1059:   for (i = 0; i < (*fvm)->numComponents; i++) {
1060:     PetscFree((*fvm)->componentNames[i]);
1061:   }
1062:   PetscFree((*fvm)->componentNames);
1063:   PetscLimiterDestroy(&(*fvm)->limiter);
1064:   PetscDualSpaceDestroy(&(*fvm)->dualSpace);
1065:   PetscFree((*fvm)->fluxWork);
1066:   PetscQuadratureDestroy(&(*fvm)->quadrature);
1067:   PetscTabulationDestroy(&(*fvm)->T);

1069:   if ((*fvm)->ops->destroy) (*(*fvm)->ops->destroy)(*fvm);
1070:   PetscHeaderDestroy(fvm);
1071:   return 0;
1072: }

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

1077:   Collective

1079:   Input Parameter:
1080: . comm - The communicator for the PetscFV object

1082:   Output Parameter:
1083: . fvm - The PetscFV object

1085:   Level: beginner

1087: .seealso: PetscFVSetType(), PETSCFVUPWIND
1088: @*/
1089: PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1090: {
1091:   PetscFV        f;

1094:   *fvm = NULL;
1095:   PetscFVInitializePackage();

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

1100:   PetscLimiterCreate(comm, &f->limiter);
1101:   f->numComponents    = 1;
1102:   f->dim              = 0;
1103:   f->computeGradients = PETSC_FALSE;
1104:   f->fluxWork         = NULL;
1105:   PetscCalloc1(f->numComponents,&f->componentNames);

1107:   *fvm = f;
1108:   return 0;
1109: }

1111: /*@
1112:   PetscFVSetLimiter - Set the limiter object

1114:   Logically collective on fvm

1116:   Input Parameters:
1117: + fvm - the PetscFV object
1118: - lim - The PetscLimiter

1120:   Level: intermediate

1122: .seealso: PetscFVGetLimiter()
1123: @*/
1124: PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1125: {
1128:   PetscLimiterDestroy(&fvm->limiter);
1129:   PetscObjectReference((PetscObject) lim);
1130:   fvm->limiter = lim;
1131:   return 0;
1132: }

1134: /*@
1135:   PetscFVGetLimiter - Get the limiter object

1137:   Not collective

1139:   Input Parameter:
1140: . fvm - the PetscFV object

1142:   Output Parameter:
1143: . lim - The PetscLimiter

1145:   Level: intermediate

1147: .seealso: PetscFVSetLimiter()
1148: @*/
1149: PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1150: {
1153:   *lim = fvm->limiter;
1154:   return 0;
1155: }

1157: /*@
1158:   PetscFVSetNumComponents - Set the number of field components

1160:   Logically collective on fvm

1162:   Input Parameters:
1163: + fvm - the PetscFV object
1164: - comp - The number of components

1166:   Level: intermediate

1168: .seealso: PetscFVGetNumComponents()
1169: @*/
1170: PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1171: {
1173:   if (fvm->numComponents != comp) {
1174:     PetscInt i;

1176:     for (i = 0; i < fvm->numComponents; i++) {
1177:       PetscFree(fvm->componentNames[i]);
1178:     }
1179:     PetscFree(fvm->componentNames);
1180:     PetscCalloc1(comp,&fvm->componentNames);
1181:   }
1182:   fvm->numComponents = comp;
1183:   PetscFree(fvm->fluxWork);
1184:   PetscMalloc1(comp, &fvm->fluxWork);
1185:   return 0;
1186: }

1188: /*@
1189:   PetscFVGetNumComponents - Get the number of field components

1191:   Not collective

1193:   Input Parameter:
1194: . fvm - the PetscFV object

1196:   Output Parameter:
1197: , comp - The number of components

1199:   Level: intermediate

1201: .seealso: PetscFVSetNumComponents()
1202: @*/
1203: PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1204: {
1207:   *comp = fvm->numComponents;
1208:   return 0;
1209: }

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

1214:   Logically collective on fvm
1215:   Input Parameters:
1216: + fvm - the PetscFV object
1217: . comp - the component number
1218: - name - the component name

1220:   Level: intermediate

1222: .seealso: PetscFVGetComponentName()
1223: @*/
1224: PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1225: {
1226:   PetscFree(fvm->componentNames[comp]);
1227:   PetscStrallocpy(name,&fvm->componentNames[comp]);
1228:   return 0;
1229: }

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

1234:   Logically collective on fvm
1235:   Input Parameters:
1236: + fvm - the PetscFV object
1237: - comp - the component number

1239:   Output Parameter:
1240: . name - the component name

1242:   Level: intermediate

1244: .seealso: PetscFVSetComponentName()
1245: @*/
1246: PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1247: {
1248:   *name = fvm->componentNames[comp];
1249:   return 0;
1250: }

1252: /*@
1253:   PetscFVSetSpatialDimension - Set the spatial dimension

1255:   Logically collective on fvm

1257:   Input Parameters:
1258: + fvm - the PetscFV object
1259: - dim - The spatial dimension

1261:   Level: intermediate

1263: .seealso: PetscFVGetSpatialDimension()
1264: @*/
1265: PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1266: {
1268:   fvm->dim = dim;
1269:   return 0;
1270: }

1272: /*@
1273:   PetscFVGetSpatialDimension - Get the spatial dimension

1275:   Logically collective on fvm

1277:   Input Parameter:
1278: . fvm - the PetscFV object

1280:   Output Parameter:
1281: . dim - The spatial dimension

1283:   Level: intermediate

1285: .seealso: PetscFVSetSpatialDimension()
1286: @*/
1287: PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1288: {
1291:   *dim = fvm->dim;
1292:   return 0;
1293: }

1295: /*@
1296:   PetscFVSetComputeGradients - Toggle computation of cell gradients

1298:   Logically collective on fvm

1300:   Input Parameters:
1301: + fvm - the PetscFV object
1302: - computeGradients - Flag to compute cell gradients

1304:   Level: intermediate

1306: .seealso: PetscFVGetComputeGradients()
1307: @*/
1308: PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1309: {
1311:   fvm->computeGradients = computeGradients;
1312:   return 0;
1313: }

1315: /*@
1316:   PetscFVGetComputeGradients - Return flag for computation of cell gradients

1318:   Not collective

1320:   Input Parameter:
1321: . fvm - the PetscFV object

1323:   Output Parameter:
1324: . computeGradients - Flag to compute cell gradients

1326:   Level: intermediate

1328: .seealso: PetscFVSetComputeGradients()
1329: @*/
1330: PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1331: {
1334:   *computeGradients = fvm->computeGradients;
1335:   return 0;
1336: }

1338: /*@
1339:   PetscFVSetQuadrature - Set the quadrature object

1341:   Logically collective on fvm

1343:   Input Parameters:
1344: + fvm - the PetscFV object
1345: - q - The PetscQuadrature

1347:   Level: intermediate

1349: .seealso: PetscFVGetQuadrature()
1350: @*/
1351: PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1352: {
1354:   PetscQuadratureDestroy(&fvm->quadrature);
1355:   PetscObjectReference((PetscObject) q);
1356:   fvm->quadrature = q;
1357:   return 0;
1358: }

1360: /*@
1361:   PetscFVGetQuadrature - Get the quadrature object

1363:   Not collective

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

1368:   Output Parameter:
1369: . lim - The PetscQuadrature

1371:   Level: intermediate

1373: .seealso: PetscFVSetQuadrature()
1374: @*/
1375: PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1376: {
1379:   if (!fvm->quadrature) {
1380:     /* Create default 1-point quadrature */
1381:     PetscReal     *points, *weights;

1383:     PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);
1384:     PetscCalloc1(fvm->dim, &points);
1385:     PetscMalloc1(1, &weights);
1386:     weights[0] = 1.0;
1387:     PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);
1388:   }
1389:   *q = fvm->quadrature;
1390:   return 0;
1391: }

1393: /*@
1394:   PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product

1396:   Not collective

1398:   Input Parameter:
1399: . fvm - The PetscFV object

1401:   Output Parameter:
1402: . sp - The PetscDualSpace object

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

1406:   Level: intermediate

1408: .seealso: PetscFVCreate()
1409: @*/
1410: PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1411: {
1414:   if (!fvm->dualSpace) {
1415:     DM              K;
1416:     PetscInt        dim, Nc, c;

1418:     PetscFVGetSpatialDimension(fvm, &dim);
1419:     PetscFVGetNumComponents(fvm, &Nc);
1420:     PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);
1421:     PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);
1422:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, PETSC_FALSE), &K);
1423:     PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);
1424:     PetscDualSpaceSetDM(fvm->dualSpace, K);
1425:     DMDestroy(&K);
1426:     PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);
1427:     /* Should we be using PetscFVGetQuadrature() here? */
1428:     for (c = 0; c < Nc; ++c) {
1429:       PetscQuadrature qc;
1430:       PetscReal      *points, *weights;

1432:       PetscQuadratureCreate(PETSC_COMM_SELF, &qc);
1433:       PetscCalloc1(dim, &points);
1434:       PetscCalloc1(Nc, &weights);
1435:       weights[c] = 1.0;
1436:       PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);
1437:       PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);
1438:       PetscQuadratureDestroy(&qc);
1439:     }
1440:     PetscDualSpaceSetUp(fvm->dualSpace);
1441:   }
1442:   *sp = fvm->dualSpace;
1443:   return 0;
1444: }

1446: /*@
1447:   PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product

1449:   Not collective

1451:   Input Parameters:
1452: + fvm - The PetscFV object
1453: - sp  - The PetscDualSpace object

1455:   Level: intermediate

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

1459: .seealso: PetscFVCreate()
1460: @*/
1461: PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1462: {
1465:   PetscDualSpaceDestroy(&fvm->dualSpace);
1466:   fvm->dualSpace = sp;
1467:   PetscObjectReference((PetscObject) fvm->dualSpace);
1468:   return 0;
1469: }

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

1474:   Not collective

1476:   Input Parameter:
1477: . fvm - The PetscFV object

1479:   Output Parameter:
1480: . T - The basis function values and derivatives at quadrature points

1482:   Note:
1483: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1484: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1485: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

1487:   Level: intermediate

1489: .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1490: @*/
1491: PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1492: {
1493:   PetscInt         npoints;
1494:   const PetscReal *points;

1498:   PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);
1499:   if (!fvm->T) PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);
1500:   *T = fvm->T;
1501:   return 0;
1502: }

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

1507:   Not collective

1509:   Input Parameters:
1510: + fvm     - The PetscFV object
1511: . nrepl   - The number of replicas
1512: . npoints - The number of tabulation points in a replica
1513: . points  - The tabulation point coordinates
1514: - K       - The order of derivative to tabulate

1516:   Output Parameter:
1517: . T - The basis function values and derivative at tabulation points

1519:   Note:
1520: $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1521: $ T->T[1] = D[((p*pdim + i)*Nc + c)*dim + d] is the derivative value at point p for basis function i, component c, in direction d
1522: $ T->T[2] = H[(((p*pdim + i)*Nc + c)*dim + d)*dim + e] is the value at point p for basis function i, component c, in directions d and e

1524:   Level: intermediate

1526: .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1527: @*/
1528: PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1529: {
1530:   PetscInt         pdim = 1; /* Dimension of approximation space P */
1531:   PetscInt         cdim;     /* Spatial dimension */
1532:   PetscInt         Nc;       /* Field components */
1533:   PetscInt         k, p, d, c, e;

1535:   if (!npoints || K < 0) {
1536:     *T = NULL;
1537:     return 0;
1538:   }
1542:   PetscFVGetSpatialDimension(fvm, &cdim);
1543:   PetscFVGetNumComponents(fvm, &Nc);
1544:   PetscMalloc1(1, T);
1545:   (*T)->K    = !cdim ? 0 : K;
1546:   (*T)->Nr   = nrepl;
1547:   (*T)->Np   = npoints;
1548:   (*T)->Nb   = pdim;
1549:   (*T)->Nc   = Nc;
1550:   (*T)->cdim = cdim;
1551:   PetscMalloc1((*T)->K+1, &(*T)->T);
1552:   for (k = 0; k <= (*T)->K; ++k) {
1553:     PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);
1554:   }
1555:   if (K >= 0) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) (*T)->T[0][(p*pdim + d)*Nc + c] = 1.0;}
1556:   if (K >= 1) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim; ++e) (*T)->T[1][((p*pdim + d)*Nc + c)*cdim + e] = 0.0;}
1557:   if (K >= 2) {for (p = 0; p < nrepl*npoints; ++p) for (d = 0; d < pdim; ++d) for (c = 0; c < Nc; ++c) for (e = 0; e < cdim*cdim; ++e) (*T)->T[2][((p*pdim + d)*Nc + c)*cdim*cdim + e] = 0.0;}
1558:   return 0;
1559: }

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

1564:   Input Parameters:
1565: + fvm      - The PetscFV object
1566: . numFaces - The number of cell faces which are not constrained
1567: - dx       - The vector from the cell centroid to the neighboring cell centroid for each face

1569:   Level: advanced

1571: .seealso: PetscFVCreate()
1572: @*/
1573: PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1574: {
1576:   if (fvm->ops->computegradient) (*fvm->ops->computegradient)(fvm, numFaces, dx, grad);
1577:   return 0;
1578: }

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

1583:   Not collective

1585:   Input Parameters:
1586: + fvm          - The PetscFV object for the field being integrated
1587: . prob         - The PetscDS specifing the discretizations and continuum functions
1588: . field        - The field being integrated
1589: . Nf           - The number of faces in the chunk
1590: . fgeom        - The face geometry for each face in the chunk
1591: . neighborVol  - The volume for each pair of cells in the chunk
1592: . uL           - The state from the cell on the left
1593: - uR           - The state from the cell on the right

1595:   Output Parameters:
1596: + fluxL        - the left fluxes for each face
1597: - fluxR        - the right fluxes for each face

1599:   Level: developer

1601: .seealso: PetscFVCreate()
1602: @*/
1603: PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1604:                                            PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1605: {
1607:   if (fvm->ops->integraterhsfunction) (*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);
1608:   return 0;
1609: }

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

1616:   Input Parameter:
1617: . fv - The initial PetscFV

1619:   Output Parameter:
1620: . fvRef - The refined PetscFV

1622:   Level: advanced

1624: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1625: @*/
1626: PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1627: {
1628:   PetscDualSpace    Q, Qref;
1629:   DM                K, Kref;
1630:   PetscQuadrature   q, qref;
1631:   DMPolytopeType    ct;
1632:   DMPlexTransform   tr;
1633:   PetscReal        *v0;
1634:   PetscReal        *jac, *invjac;
1635:   PetscInt          numComp, numSubelements, s;

1637:   PetscFVGetDualSpace(fv, &Q);
1638:   PetscFVGetQuadrature(fv, &q);
1639:   PetscDualSpaceGetDM(Q, &K);
1640:   /* Create dual space */
1641:   PetscDualSpaceDuplicate(Q, &Qref);
1642:   DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);
1643:   PetscDualSpaceSetDM(Qref, Kref);
1644:   DMDestroy(&Kref);
1645:   PetscDualSpaceSetUp(Qref);
1646:   /* Create volume */
1647:   PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);
1648:   PetscFVSetDualSpace(*fvRef, Qref);
1649:   PetscFVGetNumComponents(fv,    &numComp);
1650:   PetscFVSetNumComponents(*fvRef, numComp);
1651:   PetscFVSetUp(*fvRef);
1652:   /* Create quadrature */
1653:   DMPlexGetCellType(K, 0, &ct);
1654:   DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
1655:   DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
1656:   DMPlexRefineRegularGetAffineTransforms(tr, ct, &numSubelements, &v0, &jac, &invjac);
1657:   PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
1658:   PetscDualSpaceSimpleSetDimension(Qref, numSubelements);
1659:   for (s = 0; s < numSubelements; ++s) {
1660:     PetscQuadrature  qs;
1661:     const PetscReal *points, *weights;
1662:     PetscReal       *p, *w;
1663:     PetscInt         dim, Nc, npoints, np;

1665:     PetscQuadratureCreate(PETSC_COMM_SELF, &qs);
1666:     PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
1667:     np   = npoints/numSubelements;
1668:     PetscMalloc1(np*dim,&p);
1669:     PetscMalloc1(np*Nc,&w);
1670:     PetscArraycpy(p, &points[s*np*dim], np*dim);
1671:     PetscArraycpy(w, &weights[s*np*Nc], np*Nc);
1672:     PetscQuadratureSetData(qs, dim, Nc, np, p, w);
1673:     PetscDualSpaceSimpleSetFunctional(Qref, s, qs);
1674:     PetscQuadratureDestroy(&qs);
1675:   }
1676:   PetscFVSetQuadrature(*fvRef, qref);
1677:   DMPlexTransformDestroy(&tr);
1678:   PetscQuadratureDestroy(&qref);
1679:   PetscDualSpaceDestroy(&Qref);
1680:   return 0;
1681: }

1683: static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1684: {
1685:   PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;

1687:   PetscFree(b);
1688:   return 0;
1689: }

1691: static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1692: {
1693:   PetscInt          Nc, c;
1694:   PetscViewerFormat format;

1696:   PetscFVGetNumComponents(fv, &Nc);
1697:   PetscViewerGetFormat(viewer, &format);
1698:   PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");
1699:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1700:   for (c = 0; c < Nc; c++) {
1701:     if (fv->componentNames[c]) {
1702:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1703:     }
1704:   }
1705:   return 0;
1706: }

1708: static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1709: {
1710:   PetscBool      iascii;

1714:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1715:   if (iascii) PetscFVView_Upwind_Ascii(fv, viewer);
1716:   return 0;
1717: }

1719: static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1720: {
1721:   return 0;
1722: }

1724: /*
1725:   neighborVol[f*2+0] contains the left  geom
1726:   neighborVol[f*2+1] contains the right geom
1727: */
1728: static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1729:                                                          PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1730: {
1731:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1732:   void              *rctx;
1733:   PetscScalar       *flux = fvm->fluxWork;
1734:   const PetscScalar *constants;
1735:   PetscInt           dim, numConstants, pdim, totDim, Nc, off, f, d;

1737:   PetscDSGetTotalComponents(prob, &Nc);
1738:   PetscDSGetTotalDimension(prob, &totDim);
1739:   PetscDSGetFieldOffset(prob, field, &off);
1740:   PetscDSGetRiemannSolver(prob, field, &riemann);
1741:   PetscDSGetContext(prob, field, &rctx);
1742:   PetscDSGetConstants(prob, &numConstants, &constants);
1743:   PetscFVGetSpatialDimension(fvm, &dim);
1744:   PetscFVGetNumComponents(fvm, &pdim);
1745:   for (f = 0; f < Nf; ++f) {
1746:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1747:     for (d = 0; d < pdim; ++d) {
1748:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1749:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1750:     }
1751:   }
1752:   return 0;
1753: }

1755: static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1756: {
1757:   fvm->ops->setfromoptions          = NULL;
1758:   fvm->ops->setup                   = PetscFVSetUp_Upwind;
1759:   fvm->ops->view                    = PetscFVView_Upwind;
1760:   fvm->ops->destroy                 = PetscFVDestroy_Upwind;
1761:   fvm->ops->integraterhsfunction    = PetscFVIntegrateRHSFunction_Upwind;
1762:   return 0;
1763: }

1765: /*MC
1766:   PETSCFVUPWIND = "upwind" - A PetscFV object

1768:   Level: intermediate

1770: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1771: M*/

1773: PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1774: {
1775:   PetscFV_Upwind *b;

1778:   PetscNewLog(fvm,&b);
1779:   fvm->data = b;

1781:   PetscFVInitialize_Upwind(fvm);
1782:   return 0;
1783: }

1785: #include <petscblaslapack.h>

1787: static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1788: {
1789:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;

1791:   PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);
1792:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
1793:   PetscFree(ls);
1794:   return 0;
1795: }

1797: static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1798: {
1799:   PetscInt          Nc, c;
1800:   PetscViewerFormat format;

1802:   PetscFVGetNumComponents(fv, &Nc);
1803:   PetscViewerGetFormat(viewer, &format);
1804:   PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");
1805:   PetscViewerASCIIPrintf(viewer, "  num components: %d\n", Nc);
1806:   for (c = 0; c < Nc; c++) {
1807:     if (fv->componentNames[c]) {
1808:       PetscViewerASCIIPrintf(viewer, "    component %d: %s\n", c, fv->componentNames[c]);
1809:     }
1810:   }
1811:   return 0;
1812: }

1814: static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1815: {
1816:   PetscBool      iascii;

1820:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1821:   if (iascii) PetscFVView_LeastSquares_Ascii(fv, viewer);
1822:   return 0;
1823: }

1825: static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
1826: {
1827:   return 0;
1828: }

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

1837:   if (debug) {
1838:     PetscMalloc1(m*n,&Aback);
1839:     PetscArraycpy(Aback,A,m*n);
1840:   }

1842:   PetscBLASIntCast(m,&M);
1843:   PetscBLASIntCast(n,&N);
1844:   PetscBLASIntCast(mstride,&lda);
1845:   PetscBLASIntCast(worksize,&ldwork);
1846:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1847:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1848:   PetscFPTrapPop();
1850:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1852:   /* Extract an explicit representation of Q */
1853:   Q    = Ainv;
1854:   PetscArraycpy(Q,A,mstride*n);
1855:   K    = N;                     /* full rank */
1856:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));

1859:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1860:   Alpha = 1.0;
1861:   ldb   = lda;
1862:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1863:   /* Ainv is Q, overwritten with inverse */

1865:   if (debug) {                      /* Check that pseudo-inverse worked */
1866:     PetscScalar  Beta = 0.0;
1867:     PetscBLASInt ldc;
1868:     K   = N;
1869:     ldc = N;
1870:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1871:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1872:     PetscFree(Aback);
1873:   }
1874:   return 0;
1875: }

1877: /* Overwrites A. Can handle degenerate problems and m<n. */
1878: static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1879: {
1880:   PetscBool      debug = PETSC_FALSE;
1881:   PetscScalar   *Brhs, *Aback;
1882:   PetscScalar   *tmpwork;
1883:   PetscReal      rcond;
1884: #if defined (PETSC_USE_COMPLEX)
1885:   PetscInt       rworkSize;
1886:   PetscReal     *rwork;
1887: #endif
1888:   PetscInt       i, j, maxmn;
1889:   PetscBLASInt   M, N, lda, ldb, ldwork;
1890:   PetscBLASInt   nrhs, irank, info;

1892:   if (debug) {
1893:     PetscMalloc1(m*n,&Aback);
1894:     PetscArraycpy(Aback,A,m*n);
1895:   }

1897:   /* initialize to identity */
1898:   tmpwork = work;
1899:   Brhs = Ainv;
1900:   maxmn = PetscMax(m,n);
1901:   for (j=0; j<maxmn; j++) {
1902:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1903:   }

1905:   PetscBLASIntCast(m,&M);
1906:   PetscBLASIntCast(n,&N);
1907:   PetscBLASIntCast(mstride,&lda);
1908:   PetscBLASIntCast(maxmn,&ldb);
1909:   PetscBLASIntCast(worksize,&ldwork);
1910:   rcond = -1;
1911:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1912:   nrhs  = M;
1913: #if defined(PETSC_USE_COMPLEX)
1914:   rworkSize = 5 * PetscMin(M,N);
1915:   PetscMalloc1(rworkSize,&rwork);
1916:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info));
1917:   PetscFPTrapPop();
1918:   PetscFree(rwork);
1919: #else
1920:   nrhs  = M;
1921:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info));
1922:   PetscFPTrapPop();
1923: #endif
1925:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1927:   return 0;
1928: }

1930: #if 0
1931: static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1932: {
1933:   PetscReal       grad[2] = {0, 0};
1934:   const PetscInt *faces;
1935:   PetscInt        numFaces, f;

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

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

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

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

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

1969:   Level: developer

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

1980:   if (numFaces > maxFaces) {
1982:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
1983:   }
1984:   PetscFVGetSpatialDimension(fvm, &dim);
1985:   for (f = 0; f < numFaces; ++f) {
1986:     for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
1987:   }
1988:   /* Overwrites B with garbage, returns Binv in row-major format */
1989:   if (useSVD) {
1990:     PetscInt maxmn = PetscMax(numFaces, dim);
1991:     PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1992:     /* Binv shaped in column-major, coldim=maxmn.*/
1993:     for (f = 0; f < numFaces; ++f) {
1994:       for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d + maxmn*f];
1995:     }
1996:   } else {
1997:     PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);
1998:     /* Binv shaped in row-major, rowdim=maxFaces.*/
1999:     for (f = 0; f < numFaces; ++f) {
2000:       for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces + f];
2001:     }
2002:   }
2003:   return 0;
2004: }

2006: /*
2007:   neighborVol[f*2+0] contains the left  geom
2008:   neighborVol[f*2+1] contains the right geom
2009: */
2010: static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2011:                                                                PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2012: {
2013:   void             (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2014:   void              *rctx;
2015:   PetscScalar       *flux = fvm->fluxWork;
2016:   const PetscScalar *constants;
2017:   PetscInt           dim, numConstants, pdim, Nc, totDim, off, f, d;

2019:   PetscDSGetTotalComponents(prob, &Nc);
2020:   PetscDSGetTotalDimension(prob, &totDim);
2021:   PetscDSGetFieldOffset(prob, field, &off);
2022:   PetscDSGetRiemannSolver(prob, field, &riemann);
2023:   PetscDSGetContext(prob, field, &rctx);
2024:   PetscDSGetConstants(prob, &numConstants, &constants);
2025:   PetscFVGetSpatialDimension(fvm, &dim);
2026:   PetscFVGetNumComponents(fvm, &pdim);
2027:   for (f = 0; f < Nf; ++f) {
2028:     (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2029:     for (d = 0; d < pdim; ++d) {
2030:       fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2031:       fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2032:     }
2033:   }
2034:   return 0;
2035: }

2037: static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2038: {
2039:   PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2040:   PetscInt              dim,m,n,nrhs,minmn,maxmn;

2043:   PetscFVGetSpatialDimension(fvm, &dim);
2044:   PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);
2045:   ls->maxFaces = maxFaces;
2046:   m       = ls->maxFaces;
2047:   n       = dim;
2048:   nrhs    = ls->maxFaces;
2049:   minmn   = PetscMin(m,n);
2050:   maxmn   = PetscMax(m,n);
2051:   ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2052:   PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);
2053:   return 0;
2054: }

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

2067: /*MC
2068:   PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object

2070:   Level: intermediate

2072: .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2073: M*/

2075: PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2076: {
2077:   PetscFV_LeastSquares *ls;

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: }

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

2099:   Not collective

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

2105:   Level: intermediate

2107: .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2108: @*/
2109: PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2110: {
2112:   PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));
2113:   return 0;
2114: }