[petsc-users] Linear Reconstruction of Gradients

[Shashwat Tiwari]

Hi,
I am trying to write a second order upwind finite volume scheme on
unstructured grid using DMPlex. I am performing gradient reconstruction
using the "DMPlexReconstructGradientsFVM" function as well as a
"ComputeGradient" function that I have written myself, in which, I loop
over faces and add the gradient contribution to neighboring cells of the
face (similar to what is done in "DMPlexReconstructGradients_Internal") .
The results don't seem correct when using the
"DMPlexReconstructGradientsFVM" function where as when using my own
function, the second order is achieved. To investigate, I looked at the
gradients computed by these two functions, and noticed that
"DMPlexReconstructGradientsFVM" assigns a zero value to the gradient of
cells with less than three neighbors, for example, boundary cells (which I
assume is because leastsquares requires over-constrained system of
equations and hence, atleast three neighbors), where as, I am not
considering this factor in my "ComputeGradient" function and adding the
contribution to neighboring cells for all interior and partition faces.
Other than this, both the functions seem to compute similar values for
gradients on all the other cells and I suspect, zero values for certain
cells might be causing the problem. I wanted to know if I am missing out on
some crucial step in the code which would enable the
"DMPlexReconstructGradientsFVM" function to compute gradients for cells
with lesser number of neighbors, and also if there is a way to augment the
leastsquares stencil, i.e. to use, say, vertex neighbors instead of face
neighbors to reconstruct gradients. Also, kindly let me know if my analysis
of the problem is wrong and if there might be some other issue.

I am attaching my code for your reference. You can find the use of both
"DMPlexReconstructGradientsFVM" function and "ComputeGradient" function
(just uncomment the one you want to use) inside "ComputeResidual" function.
I am also attaching the Gmsh file for the mesh I am using. You can run the
code as follows:
          mpiexec -n <NumProc> ./convect -mesh square_tri.msh

Please let me know if you need anything else.

Thanks and Regards,
Shashwat
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200512/726b987a/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: convect.c
Type: text/x-csrc
Size: 26290 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200512/726b987a/attachment-0001.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: square_tri.geo
Type: application/octet-stream
Size: 376 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200512/726b987a/attachment-0001.obj>