Extrapolating boundary condition

[ibutsky]

Adding a new external boundary type, "extrapolate", meant to help a simulation stay in hydrostatic equilibrium. For most baryon fields, this method predicts the values in the boundary with a quadratic interpolation from the 3 nearest active cells. The density field is fit with a power law using the nearest 2 active cells. All velocities are set to zero in the boundary.

[gregbryan]

Thanks Iryna! This looks great. A few minor comments:

1. I wonder if a better name would be "hydrostatic" rather than "extrapolate"? Since you do some very specific things (use quadratic for some variables, log for others, set velocities to zero) that seem like they are really for maintaining hydrostatic equilibrium at a boundary rather than merely extrapolating.
2. Can you add a bit more description in the parameters/index.rst file to describe what the new boundary type actually does (the description for this pull request would be fine, for example).
3. I wonder if the ordering of this type and "Undefined" should be switched? I was thinking the "Undefined" should always be the last type.
   Otherwise, this looks very good to me.

[ibutsky]

Thanks @gregbryan! I implemented the changes you suggested in this new commit.

[bwoshea]

@ibutsky could you please update your pull request using the tip of the enzo-dev repository? We updated the gold standard for tests.

[jwise77]

Everything looks good to me, except that the code should either (1) halt and give an error message when using hydrostatic BCs and OOC_BOUNDARY or (2) your additions should be duplicated in those code segments.

[dcollins4096]

This looks good to me. I can see this being very useful.

Do you have a simple test problem that you could add to the run directory that uses this boundary? It's helpful to have working examples of machinery.

[jwise77]

@bwoshea Could you review this PR? Thanks!

[ibutsky]

After working with these hydrostatic boundary conditions in my thermal instability simulations, I decided to update the algorithm. I just pushed an update in commit #561b3b17ff2c1cb449ca242e32c7ee385a7a7e9d.

The biggest differences are:

1. The boundary values are now extrapolated from the slice-averaged values near the boundary. This prevents a small local over(under)density from being extrapolated to ridiculous values.
2. Only the velocity component in the direction of the boundary is set to zero, the other two components are extrapolated from average values.
3. Similarly, only the B-field component in the direction of the boundary is set to zero. Unlike the velocity, the other two components of the b-field are simply set to their values at the boundary.

[jwise77]

Hi @ibutsky, sorry for the delay in the review. But your changes don't compile correctly anymore. It looks like there's a mismatch in the routine's arguments. Could you fix this? Thanks!

[ibutsky]

@jwise77 Should be fixed now - sorry about that!

[jwise77]

Ping @bwoshea for your review. Thanks!

[jwise77]

I'm merging after we agreed on the mailing list to only require 2 approvals for PRs. Sorry this took so long!