1 #include <petsc/private/petscfvimpl.h> /*I "petscfv.h" I*/
2 #include <petscdmplex.h>
3 #include <petscds.h>
4
5 PetscClassId PETSCLIMITER_CLASSID = 0;
6
7 PetscFunctionList PetscLimiterList = NULL;
8 PetscBool PetscLimiterRegisterAllCalled = PETSC_FALSE;
9
10 PetscBool Limitercite = PETSC_FALSE;
11 const char LimiterCitation[] = "@article{BergerAftosmisMurman2005,\n"
12 " title = {Analysis of slope limiters on irregular grids},\n"
13 " journal = {AIAA paper},\n"
14 " author = {Marsha Berger and Michael J. Aftosmis and Scott M. Murman},\n"
15 " volume = {490},\n"
16 " year = {2005}\n}\n";
17
18 /*@C
19 PetscLimiterRegister - Adds a new PetscLimiter implementation
20
21 Not Collective
22
23 Input Parameters:
24 + name - The name of a new user-defined creation routine
25 - create_func - The creation routine itself
26
27 Notes:
28 PetscLimiterRegister() may be called multiple times to add several user-defined PetscLimiters
29
30 Sample usage:
31 .vb
32 PetscLimiterRegister("my_lim", MyPetscLimiterCreate);
33 .ve
34
35 Then, your PetscLimiter type can be chosen with the procedural interface via
36 .vb
37 PetscLimiterCreate(MPI_Comm, PetscLimiter *);
38 PetscLimiterSetType(PetscLimiter, "my_lim");
39 .ve
40 or at runtime via the option
41 .vb
42 -petsclimiter_type my_lim
43 .ve
44
45 Level: advanced
46
47 .seealso: PetscLimiterRegisterAll(), PetscLimiterRegisterDestroy()
48
49 @*/
PetscLimiterRegister(const char sname[],PetscErrorCode (* function)(PetscLimiter))50 PetscErrorCode PetscLimiterRegister(const char sname[], PetscErrorCode (*function)(PetscLimiter))
51 {
52 PetscErrorCode ierr;
53
54 PetscFunctionBegin;
55 ierr = PetscFunctionListAdd(&PetscLimiterList, sname, function);CHKERRQ(ierr);
56 PetscFunctionReturn(0);
57 }
58
59 /*@C
60 PetscLimiterSetType - Builds a particular PetscLimiter
61
62 Collective on lim
63
64 Input Parameters:
65 + lim - The PetscLimiter object
66 - name - The kind of limiter
67
68 Options Database Key:
69 . -petsclimiter_type <type> - Sets the PetscLimiter type; use -help for a list of available types
70
71 Level: intermediate
72
73 .seealso: PetscLimiterGetType(), PetscLimiterCreate()
74 @*/
PetscLimiterSetType(PetscLimiter lim,PetscLimiterType name)75 PetscErrorCode PetscLimiterSetType(PetscLimiter lim, PetscLimiterType name)
76 {
77 PetscErrorCode (*r)(PetscLimiter);
78 PetscBool match;
79 PetscErrorCode ierr;
80
81 PetscFunctionBegin;
82 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
83 ierr = PetscObjectTypeCompare((PetscObject) lim, name, &match);CHKERRQ(ierr);
84 if (match) PetscFunctionReturn(0);
85
86 ierr = PetscLimiterRegisterAll();CHKERRQ(ierr);
87 ierr = PetscFunctionListFind(PetscLimiterList, name, &r);CHKERRQ(ierr);
88 if (!r) SETERRQ1(PetscObjectComm((PetscObject) lim), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscLimiter type: %s", name);
89
90 if (lim->ops->destroy) {
91 ierr = (*lim->ops->destroy)(lim);CHKERRQ(ierr);
92 lim->ops->destroy = NULL;
93 }
94 ierr = (*r)(lim);CHKERRQ(ierr);
95 ierr = PetscObjectChangeTypeName((PetscObject) lim, name);CHKERRQ(ierr);
96 PetscFunctionReturn(0);
97 }
98
99 /*@C
100 PetscLimiterGetType - Gets the PetscLimiter type name (as a string) from the object.
101
102 Not Collective
103
104 Input Parameter:
105 . lim - The PetscLimiter
106
107 Output Parameter:
108 . name - The PetscLimiter type name
109
110 Level: intermediate
111
112 .seealso: PetscLimiterSetType(), PetscLimiterCreate()
113 @*/
PetscLimiterGetType(PetscLimiter lim,PetscLimiterType * name)114 PetscErrorCode PetscLimiterGetType(PetscLimiter lim, PetscLimiterType *name)
115 {
116 PetscErrorCode ierr;
117
118 PetscFunctionBegin;
119 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
120 PetscValidPointer(name, 2);
121 ierr = PetscLimiterRegisterAll();CHKERRQ(ierr);
122 *name = ((PetscObject) lim)->type_name;
123 PetscFunctionReturn(0);
124 }
125
126 /*@C
127 PetscLimiterViewFromOptions - View from Options
128
129 Collective on PetscLimiter
130
131 Input Parameters:
132 + A - the PetscLimiter object to view
133 . obj - Optional object
134 - name - command line option
135
136 Level: intermediate
137 .seealso: PetscLimiter, PetscLimiterView, PetscObjectViewFromOptions(), PetscLimiterCreate()
138 @*/
PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])139 PetscErrorCode PetscLimiterViewFromOptions(PetscLimiter A,PetscObject obj,const char name[])
140 {
141 PetscErrorCode ierr;
142
143 PetscFunctionBegin;
144 PetscValidHeaderSpecific(A,PETSCLIMITER_CLASSID,1);
145 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr);
146 PetscFunctionReturn(0);
147 }
148
149 /*@C
150 PetscLimiterView - Views a PetscLimiter
151
152 Collective on lim
153
154 Input Parameter:
155 + lim - the PetscLimiter object to view
156 - v - the viewer
157
158 Level: beginner
159
160 .seealso: PetscLimiterDestroy()
161 @*/
PetscLimiterView(PetscLimiter lim,PetscViewer v)162 PetscErrorCode PetscLimiterView(PetscLimiter lim, PetscViewer v)
163 {
164 PetscErrorCode ierr;
165
166 PetscFunctionBegin;
167 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
168 if (!v) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) lim), &v);CHKERRQ(ierr);}
169 if (lim->ops->view) {ierr = (*lim->ops->view)(lim, v);CHKERRQ(ierr);}
170 PetscFunctionReturn(0);
171 }
172
173 /*@
174 PetscLimiterSetFromOptions - sets parameters in a PetscLimiter from the options database
175
176 Collective on lim
177
178 Input Parameter:
179 . lim - the PetscLimiter object to set options for
180
181 Level: intermediate
182
183 .seealso: PetscLimiterView()
184 @*/
PetscLimiterSetFromOptions(PetscLimiter lim)185 PetscErrorCode PetscLimiterSetFromOptions(PetscLimiter lim)
186 {
187 const char *defaultType;
188 char name[256];
189 PetscBool flg;
190 PetscErrorCode ierr;
191
192 PetscFunctionBegin;
193 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
194 if (!((PetscObject) lim)->type_name) defaultType = PETSCLIMITERSIN;
195 else defaultType = ((PetscObject) lim)->type_name;
196 ierr = PetscLimiterRegisterAll();CHKERRQ(ierr);
197
198 ierr = PetscObjectOptionsBegin((PetscObject) lim);CHKERRQ(ierr);
199 ierr = PetscOptionsFList("-petsclimiter_type", "Finite volume slope limiter", "PetscLimiterSetType", PetscLimiterList, defaultType, name, 256, &flg);CHKERRQ(ierr);
200 if (flg) {
201 ierr = PetscLimiterSetType(lim, name);CHKERRQ(ierr);
202 } else if (!((PetscObject) lim)->type_name) {
203 ierr = PetscLimiterSetType(lim, defaultType);CHKERRQ(ierr);
204 }
205 if (lim->ops->setfromoptions) {ierr = (*lim->ops->setfromoptions)(lim);CHKERRQ(ierr);}
206 /* process any options handlers added with PetscObjectAddOptionsHandler() */
207 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) lim);CHKERRQ(ierr);
208 ierr = PetscOptionsEnd();CHKERRQ(ierr);
209 ierr = PetscLimiterViewFromOptions(lim, NULL, "-petsclimiter_view");CHKERRQ(ierr);
210 PetscFunctionReturn(0);
211 }
212
213 /*@C
214 PetscLimiterSetUp - Construct data structures for the PetscLimiter
215
216 Collective on lim
217
218 Input Parameter:
219 . lim - the PetscLimiter object to setup
220
221 Level: intermediate
222
223 .seealso: PetscLimiterView(), PetscLimiterDestroy()
224 @*/
PetscLimiterSetUp(PetscLimiter lim)225 PetscErrorCode PetscLimiterSetUp(PetscLimiter lim)
226 {
227 PetscErrorCode ierr;
228
229 PetscFunctionBegin;
230 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
231 if (lim->ops->setup) {ierr = (*lim->ops->setup)(lim);CHKERRQ(ierr);}
232 PetscFunctionReturn(0);
233 }
234
235 /*@
236 PetscLimiterDestroy - Destroys a PetscLimiter object
237
238 Collective on lim
239
240 Input Parameter:
241 . lim - the PetscLimiter object to destroy
242
243 Level: beginner
244
245 .seealso: PetscLimiterView()
246 @*/
PetscLimiterDestroy(PetscLimiter * lim)247 PetscErrorCode PetscLimiterDestroy(PetscLimiter *lim)
248 {
249 PetscErrorCode ierr;
250
251 PetscFunctionBegin;
252 if (!*lim) PetscFunctionReturn(0);
253 PetscValidHeaderSpecific((*lim), PETSCLIMITER_CLASSID, 1);
254
255 if (--((PetscObject)(*lim))->refct > 0) {*lim = NULL; PetscFunctionReturn(0);}
256 ((PetscObject) (*lim))->refct = 0;
257
258 if ((*lim)->ops->destroy) {ierr = (*(*lim)->ops->destroy)(*lim);CHKERRQ(ierr);}
259 ierr = PetscHeaderDestroy(lim);CHKERRQ(ierr);
260 PetscFunctionReturn(0);
261 }
262
263 /*@
264 PetscLimiterCreate - Creates an empty PetscLimiter object. The type can then be set with PetscLimiterSetType().
265
266 Collective
267
268 Input Parameter:
269 . comm - The communicator for the PetscLimiter object
270
271 Output Parameter:
272 . lim - The PetscLimiter object
273
274 Level: beginner
275
276 .seealso: PetscLimiterSetType(), PETSCLIMITERSIN
277 @*/
PetscLimiterCreate(MPI_Comm comm,PetscLimiter * lim)278 PetscErrorCode PetscLimiterCreate(MPI_Comm comm, PetscLimiter *lim)
279 {
280 PetscLimiter l;
281 PetscErrorCode ierr;
282
283 PetscFunctionBegin;
284 PetscValidPointer(lim, 2);
285 ierr = PetscCitationsRegister(LimiterCitation,&Limitercite);CHKERRQ(ierr);
286 *lim = NULL;
287 ierr = PetscFVInitializePackage();CHKERRQ(ierr);
288
289 ierr = PetscHeaderCreate(l, PETSCLIMITER_CLASSID, "PetscLimiter", "Finite Volume Slope Limiter", "PetscLimiter", comm, PetscLimiterDestroy, PetscLimiterView);CHKERRQ(ierr);
290
291 *lim = l;
292 PetscFunctionReturn(0);
293 }
294
295 /*@
296 PetscLimiterLimit - Limit the flux
297
298 Input Parameters:
299 + lim - The PetscLimiter
300 - flim - The input field
301
302 Output Parameter:
303 . phi - The limited field
304
305 Note: Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
306 $ The classical flux-limited formulation is psi(r) where
307 $
308 $ r = (u[0] - u[-1]) / (u[1] - u[0])
309 $
310 $ The second order TVD region is bounded by
311 $
312 $ psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
313 $
314 $ where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
315 $ phi(r)(r+1)/2 in which the reconstructed interface values are
316 $
317 $ u(v) = u[0] + phi(r) (grad u)[0] v
318 $
319 $ where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
320 $
321 $ 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))
322 $
323 $ For a nicer symmetric formulation, rewrite in terms of
324 $
325 $ f = (u[0] - u[-1]) / (u[1] - u[-1])
326 $
327 $ where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
328 $
329 $ phi(r) = phi(1/r)
330 $
331 $ becomes
332 $
333 $ w(f) = w(1-f).
334 $
335 $ The limiters below implement this final form w(f). The reference methods are
336 $
337 $ w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
338
339 Level: beginner
340
341 .seealso: PetscLimiterSetType(), PetscLimiterCreate()
342 @*/
PetscLimiterLimit(PetscLimiter lim,PetscReal flim,PetscReal * phi)343 PetscErrorCode PetscLimiterLimit(PetscLimiter lim, PetscReal flim, PetscReal *phi)
344 {
345 PetscErrorCode ierr;
346
347 PetscFunctionBegin;
348 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
349 PetscValidPointer(phi, 3);
350 ierr = (*lim->ops->limit)(lim, flim, phi);CHKERRQ(ierr);
351 PetscFunctionReturn(0);
352 }
353
PetscLimiterDestroy_Sin(PetscLimiter lim)354 static PetscErrorCode PetscLimiterDestroy_Sin(PetscLimiter lim)
355 {
356 PetscLimiter_Sin *l = (PetscLimiter_Sin *) lim->data;
357 PetscErrorCode ierr;
358
359 PetscFunctionBegin;
360 ierr = PetscFree(l);CHKERRQ(ierr);
361 PetscFunctionReturn(0);
362 }
363
PetscLimiterView_Sin_Ascii(PetscLimiter lim,PetscViewer viewer)364 static PetscErrorCode PetscLimiterView_Sin_Ascii(PetscLimiter lim, PetscViewer viewer)
365 {
366 PetscViewerFormat format;
367 PetscErrorCode ierr;
368
369 PetscFunctionBegin;
370 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
371 ierr = PetscViewerASCIIPrintf(viewer, "Sin Slope Limiter:\n");CHKERRQ(ierr);
372 PetscFunctionReturn(0);
373 }
374
PetscLimiterView_Sin(PetscLimiter lim,PetscViewer viewer)375 static PetscErrorCode PetscLimiterView_Sin(PetscLimiter lim, PetscViewer viewer)
376 {
377 PetscBool iascii;
378 PetscErrorCode ierr;
379
380 PetscFunctionBegin;
381 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
382 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
383 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
384 if (iascii) {ierr = PetscLimiterView_Sin_Ascii(lim, viewer);CHKERRQ(ierr);}
385 PetscFunctionReturn(0);
386 }
387
PetscLimiterLimit_Sin(PetscLimiter lim,PetscReal f,PetscReal * phi)388 static PetscErrorCode PetscLimiterLimit_Sin(PetscLimiter lim, PetscReal f, PetscReal *phi)
389 {
390 PetscFunctionBegin;
391 *phi = PetscSinReal(PETSC_PI*PetscMax(0, PetscMin(f, 1)));
392 PetscFunctionReturn(0);
393 }
394
PetscLimiterInitialize_Sin(PetscLimiter lim)395 static PetscErrorCode PetscLimiterInitialize_Sin(PetscLimiter lim)
396 {
397 PetscFunctionBegin;
398 lim->ops->view = PetscLimiterView_Sin;
399 lim->ops->destroy = PetscLimiterDestroy_Sin;
400 lim->ops->limit = PetscLimiterLimit_Sin;
401 PetscFunctionReturn(0);
402 }
403
404 /*MC
405 PETSCLIMITERSIN = "sin" - A PetscLimiter object
406
407 Level: intermediate
408
409 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
410 M*/
411
PetscLimiterCreate_Sin(PetscLimiter lim)412 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Sin(PetscLimiter lim)
413 {
414 PetscLimiter_Sin *l;
415 PetscErrorCode ierr;
416
417 PetscFunctionBegin;
418 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
419 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
420 lim->data = l;
421
422 ierr = PetscLimiterInitialize_Sin(lim);CHKERRQ(ierr);
423 PetscFunctionReturn(0);
424 }
425
PetscLimiterDestroy_Zero(PetscLimiter lim)426 static PetscErrorCode PetscLimiterDestroy_Zero(PetscLimiter lim)
427 {
428 PetscLimiter_Zero *l = (PetscLimiter_Zero *) lim->data;
429 PetscErrorCode ierr;
430
431 PetscFunctionBegin;
432 ierr = PetscFree(l);CHKERRQ(ierr);
433 PetscFunctionReturn(0);
434 }
435
PetscLimiterView_Zero_Ascii(PetscLimiter lim,PetscViewer viewer)436 static PetscErrorCode PetscLimiterView_Zero_Ascii(PetscLimiter lim, PetscViewer viewer)
437 {
438 PetscViewerFormat format;
439 PetscErrorCode ierr;
440
441 PetscFunctionBegin;
442 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
443 ierr = PetscViewerASCIIPrintf(viewer, "Zero Slope Limiter:\n");CHKERRQ(ierr);
444 PetscFunctionReturn(0);
445 }
446
PetscLimiterView_Zero(PetscLimiter lim,PetscViewer viewer)447 static PetscErrorCode PetscLimiterView_Zero(PetscLimiter lim, PetscViewer viewer)
448 {
449 PetscBool iascii;
450 PetscErrorCode ierr;
451
452 PetscFunctionBegin;
453 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
454 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
455 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
456 if (iascii) {ierr = PetscLimiterView_Zero_Ascii(lim, viewer);CHKERRQ(ierr);}
457 PetscFunctionReturn(0);
458 }
459
PetscLimiterLimit_Zero(PetscLimiter lim,PetscReal f,PetscReal * phi)460 static PetscErrorCode PetscLimiterLimit_Zero(PetscLimiter lim, PetscReal f, PetscReal *phi)
461 {
462 PetscFunctionBegin;
463 *phi = 0.0;
464 PetscFunctionReturn(0);
465 }
466
PetscLimiterInitialize_Zero(PetscLimiter lim)467 static PetscErrorCode PetscLimiterInitialize_Zero(PetscLimiter lim)
468 {
469 PetscFunctionBegin;
470 lim->ops->view = PetscLimiterView_Zero;
471 lim->ops->destroy = PetscLimiterDestroy_Zero;
472 lim->ops->limit = PetscLimiterLimit_Zero;
473 PetscFunctionReturn(0);
474 }
475
476 /*MC
477 PETSCLIMITERZERO = "zero" - A PetscLimiter object
478
479 Level: intermediate
480
481 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
482 M*/
483
PetscLimiterCreate_Zero(PetscLimiter lim)484 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Zero(PetscLimiter lim)
485 {
486 PetscLimiter_Zero *l;
487 PetscErrorCode ierr;
488
489 PetscFunctionBegin;
490 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
491 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
492 lim->data = l;
493
494 ierr = PetscLimiterInitialize_Zero(lim);CHKERRQ(ierr);
495 PetscFunctionReturn(0);
496 }
497
PetscLimiterDestroy_None(PetscLimiter lim)498 static PetscErrorCode PetscLimiterDestroy_None(PetscLimiter lim)
499 {
500 PetscLimiter_None *l = (PetscLimiter_None *) lim->data;
501 PetscErrorCode ierr;
502
503 PetscFunctionBegin;
504 ierr = PetscFree(l);CHKERRQ(ierr);
505 PetscFunctionReturn(0);
506 }
507
PetscLimiterView_None_Ascii(PetscLimiter lim,PetscViewer viewer)508 static PetscErrorCode PetscLimiterView_None_Ascii(PetscLimiter lim, PetscViewer viewer)
509 {
510 PetscViewerFormat format;
511 PetscErrorCode ierr;
512
513 PetscFunctionBegin;
514 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
515 ierr = PetscViewerASCIIPrintf(viewer, "None Slope Limiter:\n");CHKERRQ(ierr);
516 PetscFunctionReturn(0);
517 }
518
PetscLimiterView_None(PetscLimiter lim,PetscViewer viewer)519 static PetscErrorCode PetscLimiterView_None(PetscLimiter lim, PetscViewer viewer)
520 {
521 PetscBool iascii;
522 PetscErrorCode ierr;
523
524 PetscFunctionBegin;
525 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
526 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
527 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
528 if (iascii) {ierr = PetscLimiterView_None_Ascii(lim, viewer);CHKERRQ(ierr);}
529 PetscFunctionReturn(0);
530 }
531
PetscLimiterLimit_None(PetscLimiter lim,PetscReal f,PetscReal * phi)532 static PetscErrorCode PetscLimiterLimit_None(PetscLimiter lim, PetscReal f, PetscReal *phi)
533 {
534 PetscFunctionBegin;
535 *phi = 1.0;
536 PetscFunctionReturn(0);
537 }
538
PetscLimiterInitialize_None(PetscLimiter lim)539 static PetscErrorCode PetscLimiterInitialize_None(PetscLimiter lim)
540 {
541 PetscFunctionBegin;
542 lim->ops->view = PetscLimiterView_None;
543 lim->ops->destroy = PetscLimiterDestroy_None;
544 lim->ops->limit = PetscLimiterLimit_None;
545 PetscFunctionReturn(0);
546 }
547
548 /*MC
549 PETSCLIMITERNONE = "none" - A PetscLimiter object
550
551 Level: intermediate
552
553 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
554 M*/
555
PetscLimiterCreate_None(PetscLimiter lim)556 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_None(PetscLimiter lim)
557 {
558 PetscLimiter_None *l;
559 PetscErrorCode ierr;
560
561 PetscFunctionBegin;
562 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
563 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
564 lim->data = l;
565
566 ierr = PetscLimiterInitialize_None(lim);CHKERRQ(ierr);
567 PetscFunctionReturn(0);
568 }
569
PetscLimiterDestroy_Minmod(PetscLimiter lim)570 static PetscErrorCode PetscLimiterDestroy_Minmod(PetscLimiter lim)
571 {
572 PetscLimiter_Minmod *l = (PetscLimiter_Minmod *) lim->data;
573 PetscErrorCode ierr;
574
575 PetscFunctionBegin;
576 ierr = PetscFree(l);CHKERRQ(ierr);
577 PetscFunctionReturn(0);
578 }
579
PetscLimiterView_Minmod_Ascii(PetscLimiter lim,PetscViewer viewer)580 static PetscErrorCode PetscLimiterView_Minmod_Ascii(PetscLimiter lim, PetscViewer viewer)
581 {
582 PetscViewerFormat format;
583 PetscErrorCode ierr;
584
585 PetscFunctionBegin;
586 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
587 ierr = PetscViewerASCIIPrintf(viewer, "Minmod Slope Limiter:\n");CHKERRQ(ierr);
588 PetscFunctionReturn(0);
589 }
590
PetscLimiterView_Minmod(PetscLimiter lim,PetscViewer viewer)591 static PetscErrorCode PetscLimiterView_Minmod(PetscLimiter lim, PetscViewer viewer)
592 {
593 PetscBool iascii;
594 PetscErrorCode ierr;
595
596 PetscFunctionBegin;
597 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
598 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
599 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
600 if (iascii) {ierr = PetscLimiterView_Minmod_Ascii(lim, viewer);CHKERRQ(ierr);}
601 PetscFunctionReturn(0);
602 }
603
PetscLimiterLimit_Minmod(PetscLimiter lim,PetscReal f,PetscReal * phi)604 static PetscErrorCode PetscLimiterLimit_Minmod(PetscLimiter lim, PetscReal f, PetscReal *phi)
605 {
606 PetscFunctionBegin;
607 *phi = 2*PetscMax(0, PetscMin(f, 1-f));
608 PetscFunctionReturn(0);
609 }
610
PetscLimiterInitialize_Minmod(PetscLimiter lim)611 static PetscErrorCode PetscLimiterInitialize_Minmod(PetscLimiter lim)
612 {
613 PetscFunctionBegin;
614 lim->ops->view = PetscLimiterView_Minmod;
615 lim->ops->destroy = PetscLimiterDestroy_Minmod;
616 lim->ops->limit = PetscLimiterLimit_Minmod;
617 PetscFunctionReturn(0);
618 }
619
620 /*MC
621 PETSCLIMITERMINMOD = "minmod" - A PetscLimiter object
622
623 Level: intermediate
624
625 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
626 M*/
627
PetscLimiterCreate_Minmod(PetscLimiter lim)628 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Minmod(PetscLimiter lim)
629 {
630 PetscLimiter_Minmod *l;
631 PetscErrorCode ierr;
632
633 PetscFunctionBegin;
634 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
635 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
636 lim->data = l;
637
638 ierr = PetscLimiterInitialize_Minmod(lim);CHKERRQ(ierr);
639 PetscFunctionReturn(0);
640 }
641
PetscLimiterDestroy_VanLeer(PetscLimiter lim)642 static PetscErrorCode PetscLimiterDestroy_VanLeer(PetscLimiter lim)
643 {
644 PetscLimiter_VanLeer *l = (PetscLimiter_VanLeer *) lim->data;
645 PetscErrorCode ierr;
646
647 PetscFunctionBegin;
648 ierr = PetscFree(l);CHKERRQ(ierr);
649 PetscFunctionReturn(0);
650 }
651
PetscLimiterView_VanLeer_Ascii(PetscLimiter lim,PetscViewer viewer)652 static PetscErrorCode PetscLimiterView_VanLeer_Ascii(PetscLimiter lim, PetscViewer viewer)
653 {
654 PetscViewerFormat format;
655 PetscErrorCode ierr;
656
657 PetscFunctionBegin;
658 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
659 ierr = PetscViewerASCIIPrintf(viewer, "Van Leer Slope Limiter:\n");CHKERRQ(ierr);
660 PetscFunctionReturn(0);
661 }
662
PetscLimiterView_VanLeer(PetscLimiter lim,PetscViewer viewer)663 static PetscErrorCode PetscLimiterView_VanLeer(PetscLimiter lim, PetscViewer viewer)
664 {
665 PetscBool iascii;
666 PetscErrorCode ierr;
667
668 PetscFunctionBegin;
669 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
670 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
671 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
672 if (iascii) {ierr = PetscLimiterView_VanLeer_Ascii(lim, viewer);CHKERRQ(ierr);}
673 PetscFunctionReturn(0);
674 }
675
PetscLimiterLimit_VanLeer(PetscLimiter lim,PetscReal f,PetscReal * phi)676 static PetscErrorCode PetscLimiterLimit_VanLeer(PetscLimiter lim, PetscReal f, PetscReal *phi)
677 {
678 PetscFunctionBegin;
679 *phi = PetscMax(0, 4*f*(1-f));
680 PetscFunctionReturn(0);
681 }
682
PetscLimiterInitialize_VanLeer(PetscLimiter lim)683 static PetscErrorCode PetscLimiterInitialize_VanLeer(PetscLimiter lim)
684 {
685 PetscFunctionBegin;
686 lim->ops->view = PetscLimiterView_VanLeer;
687 lim->ops->destroy = PetscLimiterDestroy_VanLeer;
688 lim->ops->limit = PetscLimiterLimit_VanLeer;
689 PetscFunctionReturn(0);
690 }
691
692 /*MC
693 PETSCLIMITERVANLEER = "vanleer" - A PetscLimiter object
694
695 Level: intermediate
696
697 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
698 M*/
699
PetscLimiterCreate_VanLeer(PetscLimiter lim)700 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanLeer(PetscLimiter lim)
701 {
702 PetscLimiter_VanLeer *l;
703 PetscErrorCode ierr;
704
705 PetscFunctionBegin;
706 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
707 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
708 lim->data = l;
709
710 ierr = PetscLimiterInitialize_VanLeer(lim);CHKERRQ(ierr);
711 PetscFunctionReturn(0);
712 }
713
PetscLimiterDestroy_VanAlbada(PetscLimiter lim)714 static PetscErrorCode PetscLimiterDestroy_VanAlbada(PetscLimiter lim)
715 {
716 PetscLimiter_VanAlbada *l = (PetscLimiter_VanAlbada *) lim->data;
717 PetscErrorCode ierr;
718
719 PetscFunctionBegin;
720 ierr = PetscFree(l);CHKERRQ(ierr);
721 PetscFunctionReturn(0);
722 }
723
PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim,PetscViewer viewer)724 static PetscErrorCode PetscLimiterView_VanAlbada_Ascii(PetscLimiter lim, PetscViewer viewer)
725 {
726 PetscViewerFormat format;
727 PetscErrorCode ierr;
728
729 PetscFunctionBegin;
730 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
731 ierr = PetscViewerASCIIPrintf(viewer, "Van Albada Slope Limiter:\n");CHKERRQ(ierr);
732 PetscFunctionReturn(0);
733 }
734
PetscLimiterView_VanAlbada(PetscLimiter lim,PetscViewer viewer)735 static PetscErrorCode PetscLimiterView_VanAlbada(PetscLimiter lim, PetscViewer viewer)
736 {
737 PetscBool iascii;
738 PetscErrorCode ierr;
739
740 PetscFunctionBegin;
741 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
742 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
743 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
744 if (iascii) {ierr = PetscLimiterView_VanAlbada_Ascii(lim, viewer);CHKERRQ(ierr);}
745 PetscFunctionReturn(0);
746 }
747
PetscLimiterLimit_VanAlbada(PetscLimiter lim,PetscReal f,PetscReal * phi)748 static PetscErrorCode PetscLimiterLimit_VanAlbada(PetscLimiter lim, PetscReal f, PetscReal *phi)
749 {
750 PetscFunctionBegin;
751 *phi = PetscMax(0, 2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f)));
752 PetscFunctionReturn(0);
753 }
754
PetscLimiterInitialize_VanAlbada(PetscLimiter lim)755 static PetscErrorCode PetscLimiterInitialize_VanAlbada(PetscLimiter lim)
756 {
757 PetscFunctionBegin;
758 lim->ops->view = PetscLimiterView_VanAlbada;
759 lim->ops->destroy = PetscLimiterDestroy_VanAlbada;
760 lim->ops->limit = PetscLimiterLimit_VanAlbada;
761 PetscFunctionReturn(0);
762 }
763
764 /*MC
765 PETSCLIMITERVANALBADA = "vanalbada" - A PetscLimiter object
766
767 Level: intermediate
768
769 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
770 M*/
771
PetscLimiterCreate_VanAlbada(PetscLimiter lim)772 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_VanAlbada(PetscLimiter lim)
773 {
774 PetscLimiter_VanAlbada *l;
775 PetscErrorCode ierr;
776
777 PetscFunctionBegin;
778 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
779 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
780 lim->data = l;
781
782 ierr = PetscLimiterInitialize_VanAlbada(lim);CHKERRQ(ierr);
783 PetscFunctionReturn(0);
784 }
785
PetscLimiterDestroy_Superbee(PetscLimiter lim)786 static PetscErrorCode PetscLimiterDestroy_Superbee(PetscLimiter lim)
787 {
788 PetscLimiter_Superbee *l = (PetscLimiter_Superbee *) lim->data;
789 PetscErrorCode ierr;
790
791 PetscFunctionBegin;
792 ierr = PetscFree(l);CHKERRQ(ierr);
793 PetscFunctionReturn(0);
794 }
795
PetscLimiterView_Superbee_Ascii(PetscLimiter lim,PetscViewer viewer)796 static PetscErrorCode PetscLimiterView_Superbee_Ascii(PetscLimiter lim, PetscViewer viewer)
797 {
798 PetscViewerFormat format;
799 PetscErrorCode ierr;
800
801 PetscFunctionBegin;
802 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
803 ierr = PetscViewerASCIIPrintf(viewer, "Superbee Slope Limiter:\n");CHKERRQ(ierr);
804 PetscFunctionReturn(0);
805 }
806
PetscLimiterView_Superbee(PetscLimiter lim,PetscViewer viewer)807 static PetscErrorCode PetscLimiterView_Superbee(PetscLimiter lim, PetscViewer viewer)
808 {
809 PetscBool iascii;
810 PetscErrorCode ierr;
811
812 PetscFunctionBegin;
813 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
814 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
815 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
816 if (iascii) {ierr = PetscLimiterView_Superbee_Ascii(lim, viewer);CHKERRQ(ierr);}
817 PetscFunctionReturn(0);
818 }
819
PetscLimiterLimit_Superbee(PetscLimiter lim,PetscReal f,PetscReal * phi)820 static PetscErrorCode PetscLimiterLimit_Superbee(PetscLimiter lim, PetscReal f, PetscReal *phi)
821 {
822 PetscFunctionBegin;
823 *phi = 4*PetscMax(0, PetscMin(f, 1-f));
824 PetscFunctionReturn(0);
825 }
826
PetscLimiterInitialize_Superbee(PetscLimiter lim)827 static PetscErrorCode PetscLimiterInitialize_Superbee(PetscLimiter lim)
828 {
829 PetscFunctionBegin;
830 lim->ops->view = PetscLimiterView_Superbee;
831 lim->ops->destroy = PetscLimiterDestroy_Superbee;
832 lim->ops->limit = PetscLimiterLimit_Superbee;
833 PetscFunctionReturn(0);
834 }
835
836 /*MC
837 PETSCLIMITERSUPERBEE = "superbee" - A PetscLimiter object
838
839 Level: intermediate
840
841 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
842 M*/
843
PetscLimiterCreate_Superbee(PetscLimiter lim)844 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_Superbee(PetscLimiter lim)
845 {
846 PetscLimiter_Superbee *l;
847 PetscErrorCode ierr;
848
849 PetscFunctionBegin;
850 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
851 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
852 lim->data = l;
853
854 ierr = PetscLimiterInitialize_Superbee(lim);CHKERRQ(ierr);
855 PetscFunctionReturn(0);
856 }
857
PetscLimiterDestroy_MC(PetscLimiter lim)858 static PetscErrorCode PetscLimiterDestroy_MC(PetscLimiter lim)
859 {
860 PetscLimiter_MC *l = (PetscLimiter_MC *) lim->data;
861 PetscErrorCode ierr;
862
863 PetscFunctionBegin;
864 ierr = PetscFree(l);CHKERRQ(ierr);
865 PetscFunctionReturn(0);
866 }
867
PetscLimiterView_MC_Ascii(PetscLimiter lim,PetscViewer viewer)868 static PetscErrorCode PetscLimiterView_MC_Ascii(PetscLimiter lim, PetscViewer viewer)
869 {
870 PetscViewerFormat format;
871 PetscErrorCode ierr;
872
873 PetscFunctionBegin;
874 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
875 ierr = PetscViewerASCIIPrintf(viewer, "MC Slope Limiter:\n");CHKERRQ(ierr);
876 PetscFunctionReturn(0);
877 }
878
PetscLimiterView_MC(PetscLimiter lim,PetscViewer viewer)879 static PetscErrorCode PetscLimiterView_MC(PetscLimiter lim, PetscViewer viewer)
880 {
881 PetscBool iascii;
882 PetscErrorCode ierr;
883
884 PetscFunctionBegin;
885 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
886 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
887 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
888 if (iascii) {ierr = PetscLimiterView_MC_Ascii(lim, viewer);CHKERRQ(ierr);}
889 PetscFunctionReturn(0);
890 }
891
892 /* aka Barth-Jespersen */
PetscLimiterLimit_MC(PetscLimiter lim,PetscReal f,PetscReal * phi)893 static PetscErrorCode PetscLimiterLimit_MC(PetscLimiter lim, PetscReal f, PetscReal *phi)
894 {
895 PetscFunctionBegin;
896 *phi = PetscMin(1, 4*PetscMax(0, PetscMin(f, 1-f)));
897 PetscFunctionReturn(0);
898 }
899
PetscLimiterInitialize_MC(PetscLimiter lim)900 static PetscErrorCode PetscLimiterInitialize_MC(PetscLimiter lim)
901 {
902 PetscFunctionBegin;
903 lim->ops->view = PetscLimiterView_MC;
904 lim->ops->destroy = PetscLimiterDestroy_MC;
905 lim->ops->limit = PetscLimiterLimit_MC;
906 PetscFunctionReturn(0);
907 }
908
909 /*MC
910 PETSCLIMITERMC = "mc" - A PetscLimiter object
911
912 Level: intermediate
913
914 .seealso: PetscLimiterType, PetscLimiterCreate(), PetscLimiterSetType()
915 M*/
916
PetscLimiterCreate_MC(PetscLimiter lim)917 PETSC_EXTERN PetscErrorCode PetscLimiterCreate_MC(PetscLimiter lim)
918 {
919 PetscLimiter_MC *l;
920 PetscErrorCode ierr;
921
922 PetscFunctionBegin;
923 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 1);
924 ierr = PetscNewLog(lim, &l);CHKERRQ(ierr);
925 lim->data = l;
926
927 ierr = PetscLimiterInitialize_MC(lim);CHKERRQ(ierr);
928 PetscFunctionReturn(0);
929 }
930
931 PetscClassId PETSCFV_CLASSID = 0;
932
933 PetscFunctionList PetscFVList = NULL;
934 PetscBool PetscFVRegisterAllCalled = PETSC_FALSE;
935
936 /*@C
937 PetscFVRegister - Adds a new PetscFV implementation
938
939 Not Collective
940
941 Input Parameters:
942 + name - The name of a new user-defined creation routine
943 - create_func - The creation routine itself
944
945 Notes:
946 PetscFVRegister() may be called multiple times to add several user-defined PetscFVs
947
948 Sample usage:
949 .vb
950 PetscFVRegister("my_fv", MyPetscFVCreate);
951 .ve
952
953 Then, your PetscFV type can be chosen with the procedural interface via
954 .vb
955 PetscFVCreate(MPI_Comm, PetscFV *);
956 PetscFVSetType(PetscFV, "my_fv");
957 .ve
958 or at runtime via the option
959 .vb
960 -petscfv_type my_fv
961 .ve
962
963 Level: advanced
964
965 .seealso: PetscFVRegisterAll(), PetscFVRegisterDestroy()
966
967 @*/
PetscFVRegister(const char sname[],PetscErrorCode (* function)(PetscFV))968 PetscErrorCode PetscFVRegister(const char sname[], PetscErrorCode (*function)(PetscFV))
969 {
970 PetscErrorCode ierr;
971
972 PetscFunctionBegin;
973 ierr = PetscFunctionListAdd(&PetscFVList, sname, function);CHKERRQ(ierr);
974 PetscFunctionReturn(0);
975 }
976
977 /*@C
978 PetscFVSetType - Builds a particular PetscFV
979
980 Collective on fvm
981
982 Input Parameters:
983 + fvm - The PetscFV object
984 - name - The kind of FVM space
985
986 Options Database Key:
987 . -petscfv_type <type> - Sets the PetscFV type; use -help for a list of available types
988
989 Level: intermediate
990
991 .seealso: PetscFVGetType(), PetscFVCreate()
992 @*/
PetscFVSetType(PetscFV fvm,PetscFVType name)993 PetscErrorCode PetscFVSetType(PetscFV fvm, PetscFVType name)
994 {
995 PetscErrorCode (*r)(PetscFV);
996 PetscBool match;
997 PetscErrorCode ierr;
998
999 PetscFunctionBegin;
1000 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1001 ierr = PetscObjectTypeCompare((PetscObject) fvm, name, &match);CHKERRQ(ierr);
1002 if (match) PetscFunctionReturn(0);
1003
1004 ierr = PetscFVRegisterAll();CHKERRQ(ierr);
1005 ierr = PetscFunctionListFind(PetscFVList, name, &r);CHKERRQ(ierr);
1006 if (!r) SETERRQ1(PetscObjectComm((PetscObject) fvm), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFV type: %s", name);
1007
1008 if (fvm->ops->destroy) {
1009 ierr = (*fvm->ops->destroy)(fvm);CHKERRQ(ierr);
1010 fvm->ops->destroy = NULL;
1011 }
1012 ierr = (*r)(fvm);CHKERRQ(ierr);
1013 ierr = PetscObjectChangeTypeName((PetscObject) fvm, name);CHKERRQ(ierr);
1014 PetscFunctionReturn(0);
1015 }
1016
1017 /*@C
1018 PetscFVGetType - Gets the PetscFV type name (as a string) from the object.
1019
1020 Not Collective
1021
1022 Input Parameter:
1023 . fvm - The PetscFV
1024
1025 Output Parameter:
1026 . name - The PetscFV type name
1027
1028 Level: intermediate
1029
1030 .seealso: PetscFVSetType(), PetscFVCreate()
1031 @*/
PetscFVGetType(PetscFV fvm,PetscFVType * name)1032 PetscErrorCode PetscFVGetType(PetscFV fvm, PetscFVType *name)
1033 {
1034 PetscErrorCode ierr;
1035
1036 PetscFunctionBegin;
1037 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1038 PetscValidPointer(name, 2);
1039 ierr = PetscFVRegisterAll();CHKERRQ(ierr);
1040 *name = ((PetscObject) fvm)->type_name;
1041 PetscFunctionReturn(0);
1042 }
1043
1044 /*@C
1045 PetscFVViewFromOptions - View from Options
1046
1047 Collective on PetscFV
1048
1049 Input Parameters:
1050 + A - the PetscFV object
1051 . obj - Optional object
1052 - name - command line option
1053
1054 Level: intermediate
1055 .seealso: PetscFV, PetscFVView, PetscObjectViewFromOptions(), PetscFVCreate()
1056 @*/
PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])1057 PetscErrorCode PetscFVViewFromOptions(PetscFV A,PetscObject obj,const char name[])
1058 {
1059 PetscErrorCode ierr;
1060
1061 PetscFunctionBegin;
1062 PetscValidHeaderSpecific(A,PETSCFV_CLASSID,1);
1063 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr);
1064 PetscFunctionReturn(0);
1065 }
1066
1067 /*@C
1068 PetscFVView - Views a PetscFV
1069
1070 Collective on fvm
1071
1072 Input Parameter:
1073 + fvm - the PetscFV object to view
1074 - v - the viewer
1075
1076 Level: beginner
1077
1078 .seealso: PetscFVDestroy()
1079 @*/
PetscFVView(PetscFV fvm,PetscViewer v)1080 PetscErrorCode PetscFVView(PetscFV fvm, PetscViewer v)
1081 {
1082 PetscErrorCode ierr;
1083
1084 PetscFunctionBegin;
1085 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1086 if (!v) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fvm), &v);CHKERRQ(ierr);}
1087 if (fvm->ops->view) {ierr = (*fvm->ops->view)(fvm, v);CHKERRQ(ierr);}
1088 PetscFunctionReturn(0);
1089 }
1090
1091 /*@
1092 PetscFVSetFromOptions - sets parameters in a PetscFV from the options database
1093
1094 Collective on fvm
1095
1096 Input Parameter:
1097 . fvm - the PetscFV object to set options for
1098
1099 Options Database Key:
1100 . -petscfv_compute_gradients <bool> - Determines whether cell gradients are calculated
1101
1102 Level: intermediate
1103
1104 .seealso: PetscFVView()
1105 @*/
PetscFVSetFromOptions(PetscFV fvm)1106 PetscErrorCode PetscFVSetFromOptions(PetscFV fvm)
1107 {
1108 const char *defaultType;
1109 char name[256];
1110 PetscBool flg;
1111 PetscErrorCode ierr;
1112
1113 PetscFunctionBegin;
1114 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1115 if (!((PetscObject) fvm)->type_name) defaultType = PETSCFVUPWIND;
1116 else defaultType = ((PetscObject) fvm)->type_name;
1117 ierr = PetscFVRegisterAll();CHKERRQ(ierr);
1118
1119 ierr = PetscObjectOptionsBegin((PetscObject) fvm);CHKERRQ(ierr);
1120 ierr = PetscOptionsFList("-petscfv_type", "Finite volume discretization", "PetscFVSetType", PetscFVList, defaultType, name, 256, &flg);CHKERRQ(ierr);
1121 if (flg) {
1122 ierr = PetscFVSetType(fvm, name);CHKERRQ(ierr);
1123 } else if (!((PetscObject) fvm)->type_name) {
1124 ierr = PetscFVSetType(fvm, defaultType);CHKERRQ(ierr);
1125
1126 }
1127 ierr = PetscOptionsBool("-petscfv_compute_gradients", "Compute cell gradients", "PetscFVSetComputeGradients", fvm->computeGradients, &fvm->computeGradients, NULL);CHKERRQ(ierr);
1128 if (fvm->ops->setfromoptions) {ierr = (*fvm->ops->setfromoptions)(fvm);CHKERRQ(ierr);}
1129 /* process any options handlers added with PetscObjectAddOptionsHandler() */
1130 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fvm);CHKERRQ(ierr);
1131 ierr = PetscLimiterSetFromOptions(fvm->limiter);CHKERRQ(ierr);
1132 ierr = PetscOptionsEnd();CHKERRQ(ierr);
1133 ierr = PetscFVViewFromOptions(fvm, NULL, "-petscfv_view");CHKERRQ(ierr);
1134 PetscFunctionReturn(0);
1135 }
1136
1137 /*@
1138 PetscFVSetUp - Construct data structures for the PetscFV
1139
1140 Collective on fvm
1141
1142 Input Parameter:
1143 . fvm - the PetscFV object to setup
1144
1145 Level: intermediate
1146
1147 .seealso: PetscFVView(), PetscFVDestroy()
1148 @*/
PetscFVSetUp(PetscFV fvm)1149 PetscErrorCode PetscFVSetUp(PetscFV fvm)
1150 {
1151 PetscErrorCode ierr;
1152
1153 PetscFunctionBegin;
1154 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1155 ierr = PetscLimiterSetUp(fvm->limiter);CHKERRQ(ierr);
1156 if (fvm->ops->setup) {ierr = (*fvm->ops->setup)(fvm);CHKERRQ(ierr);}
1157 PetscFunctionReturn(0);
1158 }
1159
1160 /*@
1161 PetscFVDestroy - Destroys a PetscFV object
1162
1163 Collective on fvm
1164
1165 Input Parameter:
1166 . fvm - the PetscFV object to destroy
1167
1168 Level: beginner
1169
1170 .seealso: PetscFVView()
1171 @*/
PetscFVDestroy(PetscFV * fvm)1172 PetscErrorCode PetscFVDestroy(PetscFV *fvm)
1173 {
1174 PetscInt i;
1175 PetscErrorCode ierr;
1176
1177 PetscFunctionBegin;
1178 if (!*fvm) PetscFunctionReturn(0);
1179 PetscValidHeaderSpecific((*fvm), PETSCFV_CLASSID, 1);
1180
1181 if (--((PetscObject)(*fvm))->refct > 0) {*fvm = NULL; PetscFunctionReturn(0);}
1182 ((PetscObject) (*fvm))->refct = 0;
1183
1184 for (i = 0; i < (*fvm)->numComponents; i++) {
1185 ierr = PetscFree((*fvm)->componentNames[i]);CHKERRQ(ierr);
1186 }
1187 ierr = PetscFree((*fvm)->componentNames);CHKERRQ(ierr);
1188 ierr = PetscLimiterDestroy(&(*fvm)->limiter);CHKERRQ(ierr);
1189 ierr = PetscDualSpaceDestroy(&(*fvm)->dualSpace);CHKERRQ(ierr);
1190 ierr = PetscFree((*fvm)->fluxWork);CHKERRQ(ierr);
1191 ierr = PetscQuadratureDestroy(&(*fvm)->quadrature);CHKERRQ(ierr);
1192 ierr = PetscTabulationDestroy(&(*fvm)->T);CHKERRQ(ierr);
1193
1194 if ((*fvm)->ops->destroy) {ierr = (*(*fvm)->ops->destroy)(*fvm);CHKERRQ(ierr);}
1195 ierr = PetscHeaderDestroy(fvm);CHKERRQ(ierr);
1196 PetscFunctionReturn(0);
1197 }
1198
1199 /*@
1200 PetscFVCreate - Creates an empty PetscFV object. The type can then be set with PetscFVSetType().
1201
1202 Collective
1203
1204 Input Parameter:
1205 . comm - The communicator for the PetscFV object
1206
1207 Output Parameter:
1208 . fvm - The PetscFV object
1209
1210 Level: beginner
1211
1212 .seealso: PetscFVSetType(), PETSCFVUPWIND
1213 @*/
PetscFVCreate(MPI_Comm comm,PetscFV * fvm)1214 PetscErrorCode PetscFVCreate(MPI_Comm comm, PetscFV *fvm)
1215 {
1216 PetscFV f;
1217 PetscErrorCode ierr;
1218
1219 PetscFunctionBegin;
1220 PetscValidPointer(fvm, 2);
1221 *fvm = NULL;
1222 ierr = PetscFVInitializePackage();CHKERRQ(ierr);
1223
1224 ierr = PetscHeaderCreate(f, PETSCFV_CLASSID, "PetscFV", "Finite Volume", "PetscFV", comm, PetscFVDestroy, PetscFVView);CHKERRQ(ierr);
1225 ierr = PetscMemzero(f->ops, sizeof(struct _PetscFVOps));CHKERRQ(ierr);
1226
1227 ierr = PetscLimiterCreate(comm, &f->limiter);CHKERRQ(ierr);
1228 f->numComponents = 1;
1229 f->dim = 0;
1230 f->computeGradients = PETSC_FALSE;
1231 f->fluxWork = NULL;
1232 ierr = PetscCalloc1(f->numComponents,&f->componentNames);CHKERRQ(ierr);
1233
1234 *fvm = f;
1235 PetscFunctionReturn(0);
1236 }
1237
1238 /*@
1239 PetscFVSetLimiter - Set the limiter object
1240
1241 Logically collective on fvm
1242
1243 Input Parameters:
1244 + fvm - the PetscFV object
1245 - lim - The PetscLimiter
1246
1247 Level: intermediate
1248
1249 .seealso: PetscFVGetLimiter()
1250 @*/
PetscFVSetLimiter(PetscFV fvm,PetscLimiter lim)1251 PetscErrorCode PetscFVSetLimiter(PetscFV fvm, PetscLimiter lim)
1252 {
1253 PetscErrorCode ierr;
1254
1255 PetscFunctionBegin;
1256 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1257 PetscValidHeaderSpecific(lim, PETSCLIMITER_CLASSID, 2);
1258 ierr = PetscLimiterDestroy(&fvm->limiter);CHKERRQ(ierr);
1259 ierr = PetscObjectReference((PetscObject) lim);CHKERRQ(ierr);
1260 fvm->limiter = lim;
1261 PetscFunctionReturn(0);
1262 }
1263
1264 /*@
1265 PetscFVGetLimiter - Get the limiter object
1266
1267 Not collective
1268
1269 Input Parameter:
1270 . fvm - the PetscFV object
1271
1272 Output Parameter:
1273 . lim - The PetscLimiter
1274
1275 Level: intermediate
1276
1277 .seealso: PetscFVSetLimiter()
1278 @*/
PetscFVGetLimiter(PetscFV fvm,PetscLimiter * lim)1279 PetscErrorCode PetscFVGetLimiter(PetscFV fvm, PetscLimiter *lim)
1280 {
1281 PetscFunctionBegin;
1282 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1283 PetscValidPointer(lim, 2);
1284 *lim = fvm->limiter;
1285 PetscFunctionReturn(0);
1286 }
1287
1288 /*@
1289 PetscFVSetNumComponents - Set the number of field components
1290
1291 Logically collective on fvm
1292
1293 Input Parameters:
1294 + fvm - the PetscFV object
1295 - comp - The number of components
1296
1297 Level: intermediate
1298
1299 .seealso: PetscFVGetNumComponents()
1300 @*/
PetscFVSetNumComponents(PetscFV fvm,PetscInt comp)1301 PetscErrorCode PetscFVSetNumComponents(PetscFV fvm, PetscInt comp)
1302 {
1303 PetscErrorCode ierr;
1304
1305 PetscFunctionBegin;
1306 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1307 if (fvm->numComponents != comp) {
1308 PetscInt i;
1309
1310 for (i = 0; i < fvm->numComponents; i++) {
1311 ierr = PetscFree(fvm->componentNames[i]);CHKERRQ(ierr);
1312 }
1313 ierr = PetscFree(fvm->componentNames);CHKERRQ(ierr);
1314 ierr = PetscCalloc1(comp,&fvm->componentNames);CHKERRQ(ierr);
1315 }
1316 fvm->numComponents = comp;
1317 ierr = PetscFree(fvm->fluxWork);CHKERRQ(ierr);
1318 ierr = PetscMalloc1(comp, &fvm->fluxWork);CHKERRQ(ierr);
1319 PetscFunctionReturn(0);
1320 }
1321
1322 /*@
1323 PetscFVGetNumComponents - Get the number of field components
1324
1325 Not collective
1326
1327 Input Parameter:
1328 . fvm - the PetscFV object
1329
1330 Output Parameter:
1331 , comp - The number of components
1332
1333 Level: intermediate
1334
1335 .seealso: PetscFVSetNumComponents()
1336 @*/
PetscFVGetNumComponents(PetscFV fvm,PetscInt * comp)1337 PetscErrorCode PetscFVGetNumComponents(PetscFV fvm, PetscInt *comp)
1338 {
1339 PetscFunctionBegin;
1340 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1341 PetscValidPointer(comp, 2);
1342 *comp = fvm->numComponents;
1343 PetscFunctionReturn(0);
1344 }
1345
1346 /*@C
1347 PetscFVSetComponentName - Set the name of a component (used in output and viewing)
1348
1349 Logically collective on fvm
1350 Input Parameters:
1351 + fvm - the PetscFV object
1352 . comp - the component number
1353 - name - the component name
1354
1355 Level: intermediate
1356
1357 .seealso: PetscFVGetComponentName()
1358 @*/
PetscFVSetComponentName(PetscFV fvm,PetscInt comp,const char * name)1359 PetscErrorCode PetscFVSetComponentName(PetscFV fvm, PetscInt comp, const char *name)
1360 {
1361 PetscErrorCode ierr;
1362
1363 PetscFunctionBegin;
1364 ierr = PetscFree(fvm->componentNames[comp]);CHKERRQ(ierr);
1365 ierr = PetscStrallocpy(name,&fvm->componentNames[comp]);CHKERRQ(ierr);
1366 PetscFunctionReturn(0);
1367 }
1368
1369 /*@C
1370 PetscFVGetComponentName - Get the name of a component (used in output and viewing)
1371
1372 Logically collective on fvm
1373 Input Parameters:
1374 + fvm - the PetscFV object
1375 - comp - the component number
1376
1377 Output Parameter:
1378 . name - the component name
1379
1380 Level: intermediate
1381
1382 .seealso: PetscFVSetComponentName()
1383 @*/
PetscFVGetComponentName(PetscFV fvm,PetscInt comp,const char ** name)1384 PetscErrorCode PetscFVGetComponentName(PetscFV fvm, PetscInt comp, const char **name)
1385 {
1386 PetscFunctionBegin;
1387 *name = fvm->componentNames[comp];
1388 PetscFunctionReturn(0);
1389 }
1390
1391 /*@
1392 PetscFVSetSpatialDimension - Set the spatial dimension
1393
1394 Logically collective on fvm
1395
1396 Input Parameters:
1397 + fvm - the PetscFV object
1398 - dim - The spatial dimension
1399
1400 Level: intermediate
1401
1402 .seealso: PetscFVGetSpatialDimension()
1403 @*/
PetscFVSetSpatialDimension(PetscFV fvm,PetscInt dim)1404 PetscErrorCode PetscFVSetSpatialDimension(PetscFV fvm, PetscInt dim)
1405 {
1406 PetscFunctionBegin;
1407 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1408 fvm->dim = dim;
1409 PetscFunctionReturn(0);
1410 }
1411
1412 /*@
1413 PetscFVGetSpatialDimension - Get the spatial dimension
1414
1415 Logically collective on fvm
1416
1417 Input Parameter:
1418 . fvm - the PetscFV object
1419
1420 Output Parameter:
1421 . dim - The spatial dimension
1422
1423 Level: intermediate
1424
1425 .seealso: PetscFVSetSpatialDimension()
1426 @*/
PetscFVGetSpatialDimension(PetscFV fvm,PetscInt * dim)1427 PetscErrorCode PetscFVGetSpatialDimension(PetscFV fvm, PetscInt *dim)
1428 {
1429 PetscFunctionBegin;
1430 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1431 PetscValidPointer(dim, 2);
1432 *dim = fvm->dim;
1433 PetscFunctionReturn(0);
1434 }
1435
1436 /*@
1437 PetscFVSetComputeGradients - Toggle computation of cell gradients
1438
1439 Logically collective on fvm
1440
1441 Input Parameters:
1442 + fvm - the PetscFV object
1443 - computeGradients - Flag to compute cell gradients
1444
1445 Level: intermediate
1446
1447 .seealso: PetscFVGetComputeGradients()
1448 @*/
PetscFVSetComputeGradients(PetscFV fvm,PetscBool computeGradients)1449 PetscErrorCode PetscFVSetComputeGradients(PetscFV fvm, PetscBool computeGradients)
1450 {
1451 PetscFunctionBegin;
1452 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1453 fvm->computeGradients = computeGradients;
1454 PetscFunctionReturn(0);
1455 }
1456
1457 /*@
1458 PetscFVGetComputeGradients - Return flag for computation of cell gradients
1459
1460 Not collective
1461
1462 Input Parameter:
1463 . fvm - the PetscFV object
1464
1465 Output Parameter:
1466 . computeGradients - Flag to compute cell gradients
1467
1468 Level: intermediate
1469
1470 .seealso: PetscFVSetComputeGradients()
1471 @*/
PetscFVGetComputeGradients(PetscFV fvm,PetscBool * computeGradients)1472 PetscErrorCode PetscFVGetComputeGradients(PetscFV fvm, PetscBool *computeGradients)
1473 {
1474 PetscFunctionBegin;
1475 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1476 PetscValidPointer(computeGradients, 2);
1477 *computeGradients = fvm->computeGradients;
1478 PetscFunctionReturn(0);
1479 }
1480
1481 /*@
1482 PetscFVSetQuadrature - Set the quadrature object
1483
1484 Logically collective on fvm
1485
1486 Input Parameters:
1487 + fvm - the PetscFV object
1488 - q - The PetscQuadrature
1489
1490 Level: intermediate
1491
1492 .seealso: PetscFVGetQuadrature()
1493 @*/
PetscFVSetQuadrature(PetscFV fvm,PetscQuadrature q)1494 PetscErrorCode PetscFVSetQuadrature(PetscFV fvm, PetscQuadrature q)
1495 {
1496 PetscErrorCode ierr;
1497
1498 PetscFunctionBegin;
1499 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1500 ierr = PetscQuadratureDestroy(&fvm->quadrature);CHKERRQ(ierr);
1501 ierr = PetscObjectReference((PetscObject) q);CHKERRQ(ierr);
1502 fvm->quadrature = q;
1503 PetscFunctionReturn(0);
1504 }
1505
1506 /*@
1507 PetscFVGetQuadrature - Get the quadrature object
1508
1509 Not collective
1510
1511 Input Parameter:
1512 . fvm - the PetscFV object
1513
1514 Output Parameter:
1515 . lim - The PetscQuadrature
1516
1517 Level: intermediate
1518
1519 .seealso: PetscFVSetQuadrature()
1520 @*/
PetscFVGetQuadrature(PetscFV fvm,PetscQuadrature * q)1521 PetscErrorCode PetscFVGetQuadrature(PetscFV fvm, PetscQuadrature *q)
1522 {
1523 PetscFunctionBegin;
1524 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1525 PetscValidPointer(q, 2);
1526 if (!fvm->quadrature) {
1527 /* Create default 1-point quadrature */
1528 PetscReal *points, *weights;
1529 PetscErrorCode ierr;
1530
1531 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &fvm->quadrature);CHKERRQ(ierr);
1532 ierr = PetscCalloc1(fvm->dim, &points);CHKERRQ(ierr);
1533 ierr = PetscMalloc1(1, &weights);CHKERRQ(ierr);
1534 weights[0] = 1.0;
1535 ierr = PetscQuadratureSetData(fvm->quadrature, fvm->dim, 1, 1, points, weights);CHKERRQ(ierr);
1536 }
1537 *q = fvm->quadrature;
1538 PetscFunctionReturn(0);
1539 }
1540
1541 /*@
1542 PetscFVGetDualSpace - Returns the PetscDualSpace used to define the inner product
1543
1544 Not collective
1545
1546 Input Parameter:
1547 . fvm - The PetscFV object
1548
1549 Output Parameter:
1550 . sp - The PetscDualSpace object
1551
1552 Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1553
1554 Level: intermediate
1555
1556 .seealso: PetscFVCreate()
1557 @*/
PetscFVGetDualSpace(PetscFV fvm,PetscDualSpace * sp)1558 PetscErrorCode PetscFVGetDualSpace(PetscFV fvm, PetscDualSpace *sp)
1559 {
1560 PetscFunctionBegin;
1561 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1562 PetscValidPointer(sp, 2);
1563 if (!fvm->dualSpace) {
1564 DM K;
1565 PetscInt dim, Nc, c;
1566 PetscErrorCode ierr;
1567
1568 ierr = PetscFVGetSpatialDimension(fvm, &dim);CHKERRQ(ierr);
1569 ierr = PetscFVGetNumComponents(fvm, &Nc);CHKERRQ(ierr);
1570 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) fvm), &fvm->dualSpace);CHKERRQ(ierr);
1571 ierr = PetscDualSpaceSetType(fvm->dualSpace, PETSCDUALSPACESIMPLE);CHKERRQ(ierr);
1572 ierr = PetscDualSpaceCreateReferenceCell(fvm->dualSpace, dim, PETSC_FALSE, &K);CHKERRQ(ierr); /* TODO: The reference cell type should be held by the discretization object */
1573 ierr = PetscDualSpaceSetNumComponents(fvm->dualSpace, Nc);CHKERRQ(ierr);
1574 ierr = PetscDualSpaceSetDM(fvm->dualSpace, K);CHKERRQ(ierr);
1575 ierr = DMDestroy(&K);CHKERRQ(ierr);
1576 ierr = PetscDualSpaceSimpleSetDimension(fvm->dualSpace, Nc);CHKERRQ(ierr);
1577 /* Should we be using PetscFVGetQuadrature() here? */
1578 for (c = 0; c < Nc; ++c) {
1579 PetscQuadrature qc;
1580 PetscReal *points, *weights;
1581 PetscErrorCode ierr;
1582
1583 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &qc);CHKERRQ(ierr);
1584 ierr = PetscCalloc1(dim, &points);CHKERRQ(ierr);
1585 ierr = PetscCalloc1(Nc, &weights);CHKERRQ(ierr);
1586 weights[c] = 1.0;
1587 ierr = PetscQuadratureSetData(qc, dim, Nc, 1, points, weights);CHKERRQ(ierr);
1588 ierr = PetscDualSpaceSimpleSetFunctional(fvm->dualSpace, c, qc);CHKERRQ(ierr);
1589 ierr = PetscQuadratureDestroy(&qc);CHKERRQ(ierr);
1590 }
1591 ierr = PetscDualSpaceSetUp(fvm->dualSpace);CHKERRQ(ierr);
1592 }
1593 *sp = fvm->dualSpace;
1594 PetscFunctionReturn(0);
1595 }
1596
1597 /*@
1598 PetscFVSetDualSpace - Sets the PetscDualSpace used to define the inner product
1599
1600 Not collective
1601
1602 Input Parameters:
1603 + fvm - The PetscFV object
1604 - sp - The PetscDualSpace object
1605
1606 Level: intermediate
1607
1608 Note: A simple dual space is provided automatically, and the user typically will not need to override it.
1609
1610 .seealso: PetscFVCreate()
1611 @*/
PetscFVSetDualSpace(PetscFV fvm,PetscDualSpace sp)1612 PetscErrorCode PetscFVSetDualSpace(PetscFV fvm, PetscDualSpace sp)
1613 {
1614 PetscErrorCode ierr;
1615
1616 PetscFunctionBegin;
1617 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1618 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 2);
1619 ierr = PetscDualSpaceDestroy(&fvm->dualSpace);CHKERRQ(ierr);
1620 fvm->dualSpace = sp;
1621 ierr = PetscObjectReference((PetscObject) fvm->dualSpace);CHKERRQ(ierr);
1622 PetscFunctionReturn(0);
1623 }
1624
1625 /*@C
1626 PetscFVGetCellTabulation - Returns the tabulation of the basis functions at the quadrature points
1627
1628 Not collective
1629
1630 Input Parameter:
1631 . fvm - The PetscFV object
1632
1633 Output Parameter:
1634 . T - The basis function values and derviatives at quadrature points
1635
1636 Note:
1637 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1638 $ 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
1639 $ 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
1640
1641 Level: intermediate
1642
1643 .seealso: PetscFEGetCellTabulation(), PetscFVCreateTabulation(), PetscFVGetQuadrature(), PetscQuadratureGetData()
1644 @*/
PetscFVGetCellTabulation(PetscFV fvm,PetscTabulation * T)1645 PetscErrorCode PetscFVGetCellTabulation(PetscFV fvm, PetscTabulation *T)
1646 {
1647 PetscInt npoints;
1648 const PetscReal *points;
1649 PetscErrorCode ierr;
1650
1651 PetscFunctionBegin;
1652 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1653 PetscValidPointer(T, 2);
1654 ierr = PetscQuadratureGetData(fvm->quadrature, NULL, NULL, &npoints, &points, NULL);CHKERRQ(ierr);
1655 if (!fvm->T) {ierr = PetscFVCreateTabulation(fvm, 1, npoints, points, 1, &fvm->T);CHKERRQ(ierr);}
1656 *T = fvm->T;
1657 PetscFunctionReturn(0);
1658 }
1659
1660 /*@C
1661 PetscFVCreateTabulation - Tabulates the basis functions, and perhaps derivatives, at the points provided.
1662
1663 Not collective
1664
1665 Input Parameters:
1666 + fvm - The PetscFV object
1667 . nrepl - The number of replicas
1668 . npoints - The number of tabulation points in a replica
1669 . points - The tabulation point coordinates
1670 - K - The order of derivative to tabulate
1671
1672 Output Parameter:
1673 . T - The basis function values and derviative at tabulation points
1674
1675 Note:
1676 $ T->T[0] = B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c
1677 $ 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
1678 $ 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
1679
1680 Level: intermediate
1681
1682 .seealso: PetscFECreateTabulation(), PetscTabulationDestroy(), PetscFEGetCellTabulation()
1683 @*/
PetscFVCreateTabulation(PetscFV fvm,PetscInt nrepl,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation * T)1684 PetscErrorCode PetscFVCreateTabulation(PetscFV fvm, PetscInt nrepl, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation *T)
1685 {
1686 PetscInt pdim = 1; /* Dimension of approximation space P */
1687 PetscInt cdim; /* Spatial dimension */
1688 PetscInt Nc; /* Field components */
1689 PetscInt k, p, d, c, e;
1690 PetscErrorCode ierr;
1691
1692 PetscFunctionBegin;
1693 if (!npoints || K < 0) {
1694 *T = NULL;
1695 PetscFunctionReturn(0);
1696 }
1697 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1698 PetscValidPointer(points, 4);
1699 PetscValidPointer(T, 6);
1700 ierr = PetscFVGetSpatialDimension(fvm, &cdim);CHKERRQ(ierr);
1701 ierr = PetscFVGetNumComponents(fvm, &Nc);CHKERRQ(ierr);
1702 ierr = PetscMalloc1(1, T);CHKERRQ(ierr);
1703 (*T)->K = !cdim ? 0 : K;
1704 (*T)->Nr = nrepl;
1705 (*T)->Np = npoints;
1706 (*T)->Nb = pdim;
1707 (*T)->Nc = Nc;
1708 (*T)->cdim = cdim;
1709 ierr = PetscMalloc1((*T)->K+1, &(*T)->T);CHKERRQ(ierr);
1710 for (k = 0; k <= (*T)->K; ++k) {
1711 ierr = PetscMalloc1(nrepl*npoints*pdim*Nc*PetscPowInt(cdim, k), &(*T)->T[k]);CHKERRQ(ierr);
1712 }
1713 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;}
1714 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;}
1715 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;}
1716 PetscFunctionReturn(0);
1717 }
1718
1719 /*@C
1720 PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
1721
1722 Input Parameters:
1723 + fvm - The PetscFV object
1724 . numFaces - The number of cell faces which are not constrained
1725 - dx - The vector from the cell centroid to the neighboring cell centroid for each face
1726
1727 Level: advanced
1728
1729 .seealso: PetscFVCreate()
1730 @*/
PetscFVComputeGradient(PetscFV fvm,PetscInt numFaces,PetscScalar dx[],PetscScalar grad[])1731 PetscErrorCode PetscFVComputeGradient(PetscFV fvm, PetscInt numFaces, PetscScalar dx[], PetscScalar grad[])
1732 {
1733 PetscErrorCode ierr;
1734
1735 PetscFunctionBegin;
1736 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1737 if (fvm->ops->computegradient) {ierr = (*fvm->ops->computegradient)(fvm, numFaces, dx, grad);CHKERRQ(ierr);}
1738 PetscFunctionReturn(0);
1739 }
1740
1741 /*@C
1742 PetscFVIntegrateRHSFunction - Produce the cell residual vector for a chunk of elements by quadrature integration
1743
1744 Not collective
1745
1746 Input Parameters:
1747 + fvm - The PetscFV object for the field being integrated
1748 . prob - The PetscDS specifing the discretizations and continuum functions
1749 . field - The field being integrated
1750 . Nf - The number of faces in the chunk
1751 . fgeom - The face geometry for each face in the chunk
1752 . neighborVol - The volume for each pair of cells in the chunk
1753 . uL - The state from the cell on the left
1754 - uR - The state from the cell on the right
1755
1756 Output Parameter:
1757 + fluxL - the left fluxes for each face
1758 - fluxR - the right fluxes for each face
1759
1760 Level: developer
1761
1762 .seealso: PetscFVCreate()
1763 @*/
PetscFVIntegrateRHSFunction(PetscFV fvm,PetscDS prob,PetscInt field,PetscInt Nf,PetscFVFaceGeom * fgeom,PetscReal * neighborVol,PetscScalar uL[],PetscScalar uR[],PetscScalar fluxL[],PetscScalar fluxR[])1764 PetscErrorCode PetscFVIntegrateRHSFunction(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1765 PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1766 {
1767 PetscErrorCode ierr;
1768
1769 PetscFunctionBegin;
1770 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1771 if (fvm->ops->integraterhsfunction) {ierr = (*fvm->ops->integraterhsfunction)(fvm, prob, field, Nf, fgeom, neighborVol, uL, uR, fluxL, fluxR);CHKERRQ(ierr);}
1772 PetscFunctionReturn(0);
1773 }
1774
1775 /*@
1776 PetscFVRefine - Create a "refined" PetscFV object that refines the reference cell into smaller copies. This is typically used
1777 to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
1778 sparsity). It is also used to create an interpolation between regularly refined meshes.
1779
1780 Input Parameter:
1781 . fv - The initial PetscFV
1782
1783 Output Parameter:
1784 . fvRef - The refined PetscFV
1785
1786 Level: advanced
1787
1788 .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1789 @*/
PetscFVRefine(PetscFV fv,PetscFV * fvRef)1790 PetscErrorCode PetscFVRefine(PetscFV fv, PetscFV *fvRef)
1791 {
1792 PetscDualSpace Q, Qref;
1793 DM K, Kref;
1794 PetscQuadrature q, qref;
1795 DMPolytopeType ct;
1796 DMPlexCellRefiner cr;
1797 PetscReal *v0;
1798 PetscReal *jac, *invjac;
1799 PetscInt numComp, numSubelements, s;
1800 PetscErrorCode ierr;
1801
1802 PetscFunctionBegin;
1803 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr);
1804 ierr = PetscFVGetQuadrature(fv, &q);CHKERRQ(ierr);
1805 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr);
1806 /* Create dual space */
1807 ierr = PetscDualSpaceDuplicate(Q, &Qref);CHKERRQ(ierr);
1808 ierr = DMRefine(K, PetscObjectComm((PetscObject) fv), &Kref);CHKERRQ(ierr);
1809 ierr = PetscDualSpaceSetDM(Qref, Kref);CHKERRQ(ierr);
1810 ierr = DMDestroy(&Kref);CHKERRQ(ierr);
1811 ierr = PetscDualSpaceSetUp(Qref);CHKERRQ(ierr);
1812 /* Create volume */
1813 ierr = PetscFVCreate(PetscObjectComm((PetscObject) fv), fvRef);CHKERRQ(ierr);
1814 ierr = PetscFVSetDualSpace(*fvRef, Qref);CHKERRQ(ierr);
1815 ierr = PetscFVGetNumComponents(fv, &numComp);CHKERRQ(ierr);
1816 ierr = PetscFVSetNumComponents(*fvRef, numComp);CHKERRQ(ierr);
1817 ierr = PetscFVSetUp(*fvRef);CHKERRQ(ierr);
1818 /* Create quadrature */
1819 ierr = DMPlexGetCellType(K, 0, &ct);CHKERRQ(ierr);
1820 ierr = DMPlexCellRefinerCreate(K, &cr);CHKERRQ(ierr);
1821 ierr = DMPlexCellRefinerGetAffineTransforms(cr, ct, &numSubelements, &v0, &jac, &invjac);CHKERRQ(ierr);
1822 ierr = PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);CHKERRQ(ierr);
1823 ierr = PetscDualSpaceSimpleSetDimension(Qref, numSubelements);CHKERRQ(ierr);
1824 for (s = 0; s < numSubelements; ++s) {
1825 PetscQuadrature qs;
1826 const PetscReal *points, *weights;
1827 PetscReal *p, *w;
1828 PetscInt dim, Nc, npoints, np;
1829
1830 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &qs);CHKERRQ(ierr);
1831 ierr = PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);CHKERRQ(ierr);
1832 np = npoints/numSubelements;
1833 ierr = PetscMalloc1(np*dim,&p);CHKERRQ(ierr);
1834 ierr = PetscMalloc1(np*Nc,&w);CHKERRQ(ierr);
1835 ierr = PetscArraycpy(p, &points[s*np*dim], np*dim);CHKERRQ(ierr);
1836 ierr = PetscArraycpy(w, &weights[s*np*Nc], np*Nc);CHKERRQ(ierr);
1837 ierr = PetscQuadratureSetData(qs, dim, Nc, np, p, w);CHKERRQ(ierr);
1838 ierr = PetscDualSpaceSimpleSetFunctional(Qref, s, qs);CHKERRQ(ierr);
1839 ierr = PetscQuadratureDestroy(&qs);CHKERRQ(ierr);
1840 }
1841 ierr = PetscFVSetQuadrature(*fvRef, qref);CHKERRQ(ierr);
1842 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr);
1843 ierr = PetscQuadratureDestroy(&qref);CHKERRQ(ierr);
1844 ierr = PetscDualSpaceDestroy(&Qref);CHKERRQ(ierr);
1845 PetscFunctionReturn(0);
1846 }
1847
PetscFVDestroy_Upwind(PetscFV fvm)1848 static PetscErrorCode PetscFVDestroy_Upwind(PetscFV fvm)
1849 {
1850 PetscFV_Upwind *b = (PetscFV_Upwind *) fvm->data;
1851 PetscErrorCode ierr;
1852
1853 PetscFunctionBegin;
1854 ierr = PetscFree(b);CHKERRQ(ierr);
1855 PetscFunctionReturn(0);
1856 }
1857
PetscFVView_Upwind_Ascii(PetscFV fv,PetscViewer viewer)1858 static PetscErrorCode PetscFVView_Upwind_Ascii(PetscFV fv, PetscViewer viewer)
1859 {
1860 PetscInt Nc, c;
1861 PetscViewerFormat format;
1862 PetscErrorCode ierr;
1863
1864 PetscFunctionBegin;
1865 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1866 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
1867 ierr = PetscViewerASCIIPrintf(viewer, "Upwind Finite Volume:\n");CHKERRQ(ierr);
1868 ierr = PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);CHKERRQ(ierr);
1869 for (c = 0; c < Nc; c++) {
1870 if (fv->componentNames[c]) {
1871 ierr = PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);CHKERRQ(ierr);
1872 }
1873 }
1874 PetscFunctionReturn(0);
1875 }
1876
PetscFVView_Upwind(PetscFV fv,PetscViewer viewer)1877 static PetscErrorCode PetscFVView_Upwind(PetscFV fv, PetscViewer viewer)
1878 {
1879 PetscBool iascii;
1880 PetscErrorCode ierr;
1881
1882 PetscFunctionBegin;
1883 PetscValidHeaderSpecific(fv, PETSCFV_CLASSID, 1);
1884 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
1885 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
1886 if (iascii) {ierr = PetscFVView_Upwind_Ascii(fv, viewer);CHKERRQ(ierr);}
1887 PetscFunctionReturn(0);
1888 }
1889
PetscFVSetUp_Upwind(PetscFV fvm)1890 static PetscErrorCode PetscFVSetUp_Upwind(PetscFV fvm)
1891 {
1892 PetscFunctionBegin;
1893 PetscFunctionReturn(0);
1894 }
1895
1896 /*
1897 neighborVol[f*2+0] contains the left geom
1898 neighborVol[f*2+1] contains the right geom
1899 */
PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm,PetscDS prob,PetscInt field,PetscInt Nf,PetscFVFaceGeom * fgeom,PetscReal * neighborVol,PetscScalar uL[],PetscScalar uR[],PetscScalar fluxL[],PetscScalar fluxR[])1900 static PetscErrorCode PetscFVIntegrateRHSFunction_Upwind(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
1901 PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
1902 {
1903 void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
1904 void *rctx;
1905 PetscScalar *flux = fvm->fluxWork;
1906 const PetscScalar *constants;
1907 PetscInt dim, numConstants, pdim, totDim, Nc, off, f, d;
1908 PetscErrorCode ierr;
1909
1910 PetscFunctionBegin;
1911 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr);
1912 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1913 ierr = PetscDSGetFieldOffset(prob, field, &off);CHKERRQ(ierr);
1914 ierr = PetscDSGetRiemannSolver(prob, field, &riemann);CHKERRQ(ierr);
1915 ierr = PetscDSGetContext(prob, field, &rctx);CHKERRQ(ierr);
1916 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr);
1917 ierr = PetscFVGetSpatialDimension(fvm, &dim);CHKERRQ(ierr);
1918 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr);
1919 for (f = 0; f < Nf; ++f) {
1920 (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
1921 for (d = 0; d < pdim; ++d) {
1922 fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
1923 fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
1924 }
1925 }
1926 PetscFunctionReturn(0);
1927 }
1928
PetscFVInitialize_Upwind(PetscFV fvm)1929 static PetscErrorCode PetscFVInitialize_Upwind(PetscFV fvm)
1930 {
1931 PetscFunctionBegin;
1932 fvm->ops->setfromoptions = NULL;
1933 fvm->ops->setup = PetscFVSetUp_Upwind;
1934 fvm->ops->view = PetscFVView_Upwind;
1935 fvm->ops->destroy = PetscFVDestroy_Upwind;
1936 fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_Upwind;
1937 PetscFunctionReturn(0);
1938 }
1939
1940 /*MC
1941 PETSCFVUPWIND = "upwind" - A PetscFV object
1942
1943 Level: intermediate
1944
1945 .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
1946 M*/
1947
PetscFVCreate_Upwind(PetscFV fvm)1948 PETSC_EXTERN PetscErrorCode PetscFVCreate_Upwind(PetscFV fvm)
1949 {
1950 PetscFV_Upwind *b;
1951 PetscErrorCode ierr;
1952
1953 PetscFunctionBegin;
1954 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
1955 ierr = PetscNewLog(fvm,&b);CHKERRQ(ierr);
1956 fvm->data = b;
1957
1958 ierr = PetscFVInitialize_Upwind(fvm);CHKERRQ(ierr);
1959 PetscFunctionReturn(0);
1960 }
1961
1962 #include <petscblaslapack.h>
1963
PetscFVDestroy_LeastSquares(PetscFV fvm)1964 static PetscErrorCode PetscFVDestroy_LeastSquares(PetscFV fvm)
1965 {
1966 PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
1967 PetscErrorCode ierr;
1968
1969 PetscFunctionBegin;
1970 ierr = PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", NULL);CHKERRQ(ierr);
1971 ierr = PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);CHKERRQ(ierr);
1972 ierr = PetscFree(ls);CHKERRQ(ierr);
1973 PetscFunctionReturn(0);
1974 }
1975
PetscFVView_LeastSquares_Ascii(PetscFV fv,PetscViewer viewer)1976 static PetscErrorCode PetscFVView_LeastSquares_Ascii(PetscFV fv, PetscViewer viewer)
1977 {
1978 PetscInt Nc, c;
1979 PetscViewerFormat format;
1980 PetscErrorCode ierr;
1981
1982 PetscFunctionBegin;
1983 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1984 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
1985 ierr = PetscViewerASCIIPrintf(viewer, "Finite Volume with Least Squares Reconstruction:\n");CHKERRQ(ierr);
1986 ierr = PetscViewerASCIIPrintf(viewer, " num components: %d\n", Nc);CHKERRQ(ierr);
1987 for (c = 0; c < Nc; c++) {
1988 if (fv->componentNames[c]) {
1989 ierr = PetscViewerASCIIPrintf(viewer, " component %d: %s\n", c, fv->componentNames[c]);CHKERRQ(ierr);
1990 }
1991 }
1992 PetscFunctionReturn(0);
1993 }
1994
PetscFVView_LeastSquares(PetscFV fv,PetscViewer viewer)1995 static PetscErrorCode PetscFVView_LeastSquares(PetscFV fv, PetscViewer viewer)
1996 {
1997 PetscBool iascii;
1998 PetscErrorCode ierr;
1999
2000 PetscFunctionBegin;
2001 PetscValidHeaderSpecific(fv, PETSCFV_CLASSID, 1);
2002 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
2003 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
2004 if (iascii) {ierr = PetscFVView_LeastSquares_Ascii(fv, viewer);CHKERRQ(ierr);}
2005 PetscFunctionReturn(0);
2006 }
2007
PetscFVSetUp_LeastSquares(PetscFV fvm)2008 static PetscErrorCode PetscFVSetUp_LeastSquares(PetscFV fvm)
2009 {
2010 PetscFunctionBegin;
2011 PetscFunctionReturn(0);
2012 }
2013
2014 /* Overwrites A. Can only handle full-rank problems with m>=n */
PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar * A,PetscScalar * Ainv,PetscScalar * tau,PetscInt worksize,PetscScalar * work)2015 static PetscErrorCode PetscFVLeastSquaresPseudoInverse_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2016 {
2017 PetscBool debug = PETSC_FALSE;
2018 PetscErrorCode ierr;
2019 PetscBLASInt M,N,K,lda,ldb,ldwork,info;
2020 PetscScalar *R,*Q,*Aback,Alpha;
2021
2022 PetscFunctionBegin;
2023 if (debug) {
2024 ierr = PetscMalloc1(m*n,&Aback);CHKERRQ(ierr);
2025 ierr = PetscArraycpy(Aback,A,m*n);CHKERRQ(ierr);
2026 }
2027
2028 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr);
2029 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr);
2030 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr);
2031 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr);
2032 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
2033 LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
2034 ierr = PetscFPTrapPop();CHKERRQ(ierr);
2035 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
2036 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
2037
2038 /* Extract an explicit representation of Q */
2039 Q = Ainv;
2040 ierr = PetscArraycpy(Q,A,mstride*n);CHKERRQ(ierr);
2041 K = N; /* full rank */
2042 PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
2043 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
2044
2045 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
2046 Alpha = 1.0;
2047 ldb = lda;
2048 BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
2049 /* Ainv is Q, overwritten with inverse */
2050
2051 if (debug) { /* Check that pseudo-inverse worked */
2052 PetscScalar Beta = 0.0;
2053 PetscBLASInt ldc;
2054 K = N;
2055 ldc = N;
2056 BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
2057 ierr = PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
2058 ierr = PetscFree(Aback);CHKERRQ(ierr);
2059 }
2060 PetscFunctionReturn(0);
2061 }
2062
2063 /* Overwrites A. Can handle degenerate problems and m<n. */
PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar * A,PetscScalar * Ainv,PetscScalar * tau,PetscInt worksize,PetscScalar * work)2064 static PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD_Static(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
2065 {
2066 PetscBool debug = PETSC_FALSE;
2067 PetscScalar *Brhs, *Aback;
2068 PetscScalar *tmpwork;
2069 PetscReal rcond;
2070 #if defined (PETSC_USE_COMPLEX)
2071 PetscInt rworkSize;
2072 PetscReal *rwork;
2073 #endif
2074 PetscInt i, j, maxmn;
2075 PetscBLASInt M, N, lda, ldb, ldwork;
2076 PetscBLASInt nrhs, irank, info;
2077 PetscErrorCode ierr;
2078
2079 PetscFunctionBegin;
2080 if (debug) {
2081 ierr = PetscMalloc1(m*n,&Aback);CHKERRQ(ierr);
2082 ierr = PetscArraycpy(Aback,A,m*n);CHKERRQ(ierr);
2083 }
2084
2085 /* initialize to identity */
2086 tmpwork = work;
2087 Brhs = Ainv;
2088 maxmn = PetscMax(m,n);
2089 for (j=0; j<maxmn; j++) {
2090 for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
2091 }
2092
2093 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr);
2094 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr);
2095 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr);
2096 ierr = PetscBLASIntCast(maxmn,&ldb);CHKERRQ(ierr);
2097 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr);
2098 rcond = -1;
2099 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
2100 nrhs = M;
2101 #if defined(PETSC_USE_COMPLEX)
2102 rworkSize = 5 * PetscMin(M,N);
2103 ierr = PetscMalloc1(rworkSize,&rwork);CHKERRQ(ierr);
2104 LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,rwork,&info);
2105 ierr = PetscFPTrapPop();CHKERRQ(ierr);
2106 ierr = PetscFree(rwork);CHKERRQ(ierr);
2107 #else
2108 nrhs = M;
2109 LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
2110 ierr = PetscFPTrapPop();CHKERRQ(ierr);
2111 #endif
2112 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
2113 /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
2114 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");
2115 PetscFunctionReturn(0);
2116 }
2117
2118 #if 0
2119 static PetscErrorCode PetscFVLeastSquaresDebugCell_Static(PetscFV fvm, PetscInt cell, DM dm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2120 {
2121 PetscReal grad[2] = {0, 0};
2122 const PetscInt *faces;
2123 PetscInt numFaces, f;
2124 PetscErrorCode ierr;
2125
2126 PetscFunctionBegin;
2127 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr);
2128 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr);
2129 for (f = 0; f < numFaces; ++f) {
2130 const PetscInt *fcells;
2131 const CellGeom *cg1;
2132 const FaceGeom *fg;
2133
2134 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr);
2135 ierr = DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr);
2136 for (i = 0; i < 2; ++i) {
2137 PetscScalar du;
2138
2139 if (fcells[i] == c) continue;
2140 ierr = DMPlexPointLocalRead(dmCell, fcells[i], cgeom, &cg1);CHKERRQ(ierr);
2141 du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
2142 grad[0] += fg->grad[!i][0] * du;
2143 grad[1] += fg->grad[!i][1] * du;
2144 }
2145 }
2146 PetscPrintf(PETSC_COMM_SELF, "cell[%d] grad (%g, %g)\n", cell, grad[0], grad[1]);CHKERRQ(ierr);
2147 PetscFunctionReturn(0);
2148 }
2149 #endif
2150
2151 /*
2152 PetscFVComputeGradient - Compute the gradient reconstruction matrix for a given cell
2153
2154 Input Parameters:
2155 + fvm - The PetscFV object
2156 . numFaces - The number of cell faces which are not constrained
2157 . dx - The vector from the cell centroid to the neighboring cell centroid for each face
2158
2159 Level: developer
2160
2161 .seealso: PetscFVCreate()
2162 */
PetscFVComputeGradient_LeastSquares(PetscFV fvm,PetscInt numFaces,const PetscScalar dx[],PetscScalar grad[])2163 static PetscErrorCode PetscFVComputeGradient_LeastSquares(PetscFV fvm, PetscInt numFaces, const PetscScalar dx[], PetscScalar grad[])
2164 {
2165 PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2166 const PetscBool useSVD = PETSC_TRUE;
2167 const PetscInt maxFaces = ls->maxFaces;
2168 PetscInt dim, f, d;
2169 PetscErrorCode ierr;
2170
2171 PetscFunctionBegin;
2172 if (numFaces > maxFaces) {
2173 if (maxFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Reconstruction has not been initialized, call PetscFVLeastSquaresSetMaxFaces()");
2174 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input faces %D > %D maxfaces", numFaces, maxFaces);
2175 }
2176 ierr = PetscFVGetSpatialDimension(fvm, &dim);CHKERRQ(ierr);
2177 for (f = 0; f < numFaces; ++f) {
2178 for (d = 0; d < dim; ++d) ls->B[d*maxFaces+f] = dx[f*dim+d];
2179 }
2180 /* Overwrites B with garbage, returns Binv in row-major format */
2181 if (useSVD) {
2182 PetscInt maxmn = PetscMax(numFaces, dim);
2183 ierr = PetscFVLeastSquaresPseudoInverseSVD_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);CHKERRQ(ierr);
2184 /* Binv shaped in column-major, coldim=maxmn.*/
2185 for (f = 0; f < numFaces; ++f) {
2186 for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d + maxmn*f];
2187 }
2188 } else {
2189 ierr = PetscFVLeastSquaresPseudoInverse_Static(numFaces, maxFaces, dim, ls->B, ls->Binv, ls->tau, ls->workSize, ls->work);CHKERRQ(ierr);
2190 /* Binv shaped in row-major, rowdim=maxFaces.*/
2191 for (f = 0; f < numFaces; ++f) {
2192 for (d = 0; d < dim; ++d) grad[f*dim+d] = ls->Binv[d*maxFaces + f];
2193 }
2194 }
2195 PetscFunctionReturn(0);
2196 }
2197
2198 /*
2199 neighborVol[f*2+0] contains the left geom
2200 neighborVol[f*2+1] contains the right geom
2201 */
PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm,PetscDS prob,PetscInt field,PetscInt Nf,PetscFVFaceGeom * fgeom,PetscReal * neighborVol,PetscScalar uL[],PetscScalar uR[],PetscScalar fluxL[],PetscScalar fluxR[])2202 static PetscErrorCode PetscFVIntegrateRHSFunction_LeastSquares(PetscFV fvm, PetscDS prob, PetscInt field, PetscInt Nf, PetscFVFaceGeom *fgeom, PetscReal *neighborVol,
2203 PetscScalar uL[], PetscScalar uR[], PetscScalar fluxL[], PetscScalar fluxR[])
2204 {
2205 void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
2206 void *rctx;
2207 PetscScalar *flux = fvm->fluxWork;
2208 const PetscScalar *constants;
2209 PetscInt dim, numConstants, pdim, Nc, totDim, off, f, d;
2210 PetscErrorCode ierr;
2211
2212 PetscFunctionBegin;
2213 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr);
2214 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
2215 ierr = PetscDSGetFieldOffset(prob, field, &off);CHKERRQ(ierr);
2216 ierr = PetscDSGetRiemannSolver(prob, field, &riemann);CHKERRQ(ierr);
2217 ierr = PetscDSGetContext(prob, field, &rctx);CHKERRQ(ierr);
2218 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr);
2219 ierr = PetscFVGetSpatialDimension(fvm, &dim);CHKERRQ(ierr);
2220 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr);
2221 for (f = 0; f < Nf; ++f) {
2222 (*riemann)(dim, pdim, fgeom[f].centroid, fgeom[f].normal, &uL[f*Nc], &uR[f*Nc], numConstants, constants, flux, rctx);
2223 for (d = 0; d < pdim; ++d) {
2224 fluxL[f*totDim+off+d] = flux[d] / neighborVol[f*2+0];
2225 fluxR[f*totDim+off+d] = flux[d] / neighborVol[f*2+1];
2226 }
2227 }
2228 PetscFunctionReturn(0);
2229 }
2230
PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm,PetscInt maxFaces)2231 static PetscErrorCode PetscFVLeastSquaresSetMaxFaces_LS(PetscFV fvm, PetscInt maxFaces)
2232 {
2233 PetscFV_LeastSquares *ls = (PetscFV_LeastSquares *) fvm->data;
2234 PetscInt dim,m,n,nrhs,minmn,maxmn;
2235 PetscErrorCode ierr;
2236
2237 PetscFunctionBegin;
2238 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
2239 ierr = PetscFVGetSpatialDimension(fvm, &dim);CHKERRQ(ierr);
2240 ierr = PetscFree4(ls->B, ls->Binv, ls->tau, ls->work);CHKERRQ(ierr);
2241 ls->maxFaces = maxFaces;
2242 m = ls->maxFaces;
2243 n = dim;
2244 nrhs = ls->maxFaces;
2245 minmn = PetscMin(m,n);
2246 maxmn = PetscMax(m,n);
2247 ls->workSize = 3*minmn + PetscMax(2*minmn, PetscMax(maxmn, nrhs)); /* required by LAPACK */
2248 ierr = PetscMalloc4(m*n,&ls->B,maxmn*maxmn,&ls->Binv,minmn,&ls->tau,ls->workSize,&ls->work);CHKERRQ(ierr);
2249 PetscFunctionReturn(0);
2250 }
2251
PetscFVInitialize_LeastSquares(PetscFV fvm)2252 PetscErrorCode PetscFVInitialize_LeastSquares(PetscFV fvm)
2253 {
2254 PetscFunctionBegin;
2255 fvm->ops->setfromoptions = NULL;
2256 fvm->ops->setup = PetscFVSetUp_LeastSquares;
2257 fvm->ops->view = PetscFVView_LeastSquares;
2258 fvm->ops->destroy = PetscFVDestroy_LeastSquares;
2259 fvm->ops->computegradient = PetscFVComputeGradient_LeastSquares;
2260 fvm->ops->integraterhsfunction = PetscFVIntegrateRHSFunction_LeastSquares;
2261 PetscFunctionReturn(0);
2262 }
2263
2264 /*MC
2265 PETSCFVLEASTSQUARES = "leastsquares" - A PetscFV object
2266
2267 Level: intermediate
2268
2269 .seealso: PetscFVType, PetscFVCreate(), PetscFVSetType()
2270 M*/
2271
PetscFVCreate_LeastSquares(PetscFV fvm)2272 PETSC_EXTERN PetscErrorCode PetscFVCreate_LeastSquares(PetscFV fvm)
2273 {
2274 PetscFV_LeastSquares *ls;
2275 PetscErrorCode ierr;
2276
2277 PetscFunctionBegin;
2278 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
2279 ierr = PetscNewLog(fvm, &ls);CHKERRQ(ierr);
2280 fvm->data = ls;
2281
2282 ls->maxFaces = -1;
2283 ls->workSize = -1;
2284 ls->B = NULL;
2285 ls->Binv = NULL;
2286 ls->tau = NULL;
2287 ls->work = NULL;
2288
2289 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr);
2290 ierr = PetscFVInitialize_LeastSquares(fvm);CHKERRQ(ierr);
2291 ierr = PetscObjectComposeFunction((PetscObject) fvm, "PetscFVLeastSquaresSetMaxFaces_C", PetscFVLeastSquaresSetMaxFaces_LS);CHKERRQ(ierr);
2292 PetscFunctionReturn(0);
2293 }
2294
2295 /*@
2296 PetscFVLeastSquaresSetMaxFaces - Set the maximum number of cell faces for gradient reconstruction
2297
2298 Not collective
2299
2300 Input parameters:
2301 + fvm - The PetscFV object
2302 - maxFaces - The maximum number of cell faces
2303
2304 Level: intermediate
2305
2306 .seealso: PetscFVCreate(), PETSCFVLEASTSQUARES
2307 @*/
PetscFVLeastSquaresSetMaxFaces(PetscFV fvm,PetscInt maxFaces)2308 PetscErrorCode PetscFVLeastSquaresSetMaxFaces(PetscFV fvm, PetscInt maxFaces)
2309 {
2310 PetscErrorCode ierr;
2311
2312 PetscFunctionBegin;
2313 PetscValidHeaderSpecific(fvm, PETSCFV_CLASSID, 1);
2314 ierr = PetscTryMethod(fvm, "PetscFVLeastSquaresSetMaxFaces_C", (PetscFV,PetscInt), (fvm,maxFaces));CHKERRQ(ierr);
2315 PetscFunctionReturn(0);
2316 }
2317