FV: add support for higher-order solution reconstruction

[rwcarlsen]

Reason

Various numerical flux methods and other FV bells-and-whistles benefit from or require having gradient approximations for mesh cells (i.e. _grad_u_left and _grad_u_right in fv kernels).

Design

Computing these gradients requires values from a cell's neighbors. This calculation should probably be initiated on-demand for FV variables when they are computing/setting the rest of the solution, gradient, etc. values. Then the variable could provide these values to kernels, bcs, etc. transparently just like objects usually get gradients. The FV variable could cash this value for the next time it is needed (since each cell's gradient info will be needed for every visit to a face of that cell). We will need to create a new "reconstruction" object that users can write custom versions of and enable/swap from input file syntax.

Impact

Should not affect existing APIs

[YaqiWang]

You may just use higher order variable to store this instead of cache. While FV may just work on the low order part and you can have a posprocessing function to recover the high order part. Just suggesting you an option. Cache is typically hard to manage and not transparent to users.

[lindsayad]

@rwcarlsen to me there are couple different concepts here. There's the ability to (re)construct a gradient and then there's also the ability to construct a higher-order (e.g. non-constant) representation of the solution within a cell. #15825 has the gradient (re)construction but it does not do anything to construct a higher-order representation of the solution within a cell. That higher-order cell representation would be really useful I think at things like outflow boundaries where I do not simply want to use the cell center value as the face value, which would be implying a potentially inaccurate zero-gradient condition.

[rwcarlsen]

Yeah - my description here really is about gradient reconstruction - even though the title is about solution reconstruction. They are related - obviously a higher order solution needs a reconstructed gradient. I'm just not sure that an explicit higher order solution reconstruction is the right angle to attack this from. With FV - everything revolves around numerical fluxes - so I would expect the numerical flux methods to incorporate the reconstructed gradient to compute whatever equivalent higher-order solution they might need. I think it might be better to have some helper utils/functions that make it easy for objects to get a "higher order" solution - utils that just use the reconstructed gradient. I think trying to do things through the variable objects would be a bit premature perhaps - at least until we have multiple concrete uses of this in the wild that could benefit from more than just some utils/funcs.