BoundaryOp refactor

[johnomotani]

Simplify boundary_standard.cxx:

• Most boundary conditions now use BoundaryOp::apply, which handles looping and staggering. The derived classes BoundaryNeumann, etc. mostly provide just methods to apply the boundary condition at a point, e.g. https://github.com/boutproject/BOUT-dev/blob/cf7aa24b0ee2870623cb398fb22dd81b34b7b0a8/src/mesh/boundary_standard.cxx#L415-L445 BoundaryDirichlet still has its own apply() method because it handled 2nd (and 3rd, 4th, etc.) guard cells in a non-standard way, setting them equal to the boundary value expression (evaluated at the guard cell's location); I have left this as-is. BoundaryNeumann_NonOrthogonal, BoundaryZeroLaplace, BoundaryZeroLaplace2 and BoundaryConstLaplace also implement overriding apply methods.
• Remove BoundaryOpBase class. This is replaced by using templates to implement several methods for both BoundaryOp and BoundaryOpPar with a minimum of code duplication. Not totally sure this form is better/clearer. It does avoid several static_casts. Could undo this, but would take a little digging because I didn't do a good job of implementing this refactor in separate commits.

Fixes a few bugs: setting boundary conditions for fields at CELL_ZLOW (previously this case fell though and the guard cells were not changed); use the fieldmesh of the field, not the global mesh, in all boundary conditions; correct one of the coefficients in BoundaryDirichlet_4thOrder (now renamed to BoundaryDirichlet_O5; implement BoundaryNeumann_O4 for staggered grids.

Add GlobalZ method to Mesh for consistency with GlobalX, GlobalY.

Add a couple of higher order extrapolating boundary conditions, free_o4 and free_o5.

Change BoundaryRobin to take account of grid spacing, so now the derivative coefficient passed in (b) does not have to include dx or dy. Effectively derivative in boundary condition is now DDX or DDY instead of indexDDX or indexDDY.

Test the order of accuracy of Dirichlet, Neumann and free boundary conditions in MMS/derivatives3.

Add new boundary condition dirichlet_smooth, which suppresses grid-scale oscillations at the boundary enough to be used in wave-1d/wave-1d-y MMS tests. It has the same accuracy as the standard dirichlet condition, but uses two grid points to set the boundary condition. Setting the guard cell g by extrapolating from the boundary value, v, and first grid cell, f1 would give 'g1 = 2v-f1'; setting from the second grid cell would give 'g2 = 4/3v-1/3f2'; this bc takes the average of the two, setting 'g = 4/3v-1/2f1-1/6f2'.

Test both dirichlet_smooth and neumann boundary conditions for StaggerGrids=true or false in MMS/wave-1d and MMS/wave-1d-y.

Test both dirichlet and neumann boundary conditions on both boundaries in MMS/diffusion.

[johnomotani]

This should replace #1319.

[dschwoerer]

Change BoundaryRobin to take account of grid spacing, so now the derivative coefficient passed in (b) does not have to include dx or dy. Effectively derivative in boundary condition is now DDX or DDY instead of indexDDX or indexDDY.

Does it take only a uniform dx/dy into account, or also a non-uniform dx/dy?

I am asking, because to my understanding the others (I had a look at free, dirichlet, neumann) do not as far as I can tell, and that is why I added #1179

[johnomotani]

Does it take only a uniform dx/dy into account, or also a non-uniform dx/dy?

@dschwoerer no, the Robin bc doesn't account for variation of dx/dy, just uses the value dx/dy at the boundary grid point. Sorry, I hadn't looked at #1179 before now.

[d7919]

Sorry for the duplicate comments!

[johnomotani]

@bendudson what is the reason for having the BoundaryWidth modifier? It's being trickier than I'd anticipated to update; understanding the purpose would help make sure I don't break anything...

[bendudson]

@johnomotani I'm a bit unsure what the purpose was of the width modifier, and since it doesn't seem to be used in any of the examples, it can probably be removed. At least I wouldn't worry too much about trying to preserve it. I think it was either for:

• Operators like C4
• Staggered cells, where in the simple implementation the boundary needs to apply to an extra point. This extra point is now handled properly (?) by the boundary condition operators.

[johnomotani]

Thanks @bendudson :+1:

The tricky part about width keyword or BoundaryWidth modifier is that the BoundaryRegion needs to have its width changed. width is used in e.g. BoundaryRegionXIn::first(), so it's not enough to just replace uses of bndry->width in boundary_standard.cxx.

The neatest solution I've thought of is to use the BoundaryWidth modifier and have it make and store a copy of the BoundaryRegion with a different width. The downside is that this required adding a copy(int wid) method to BoundaryRegion. Suggestions for improvements or better designs welcome... or maybe we should just get rid of BoundaryWidth...

A few other changes:

• I noticed FieldData::setBoundary(name) just added boundary conditions to the bndry_op vector without checking it was empty first. Probably it only ever gets called once, but just in case 1f588a4 makes it delete the entries and clear the vector. If it was called more than once, multiple boundary conditions could be applied on the same boundary. The last one would overwrite the others, so should be correct but inefficient.

• I store the BoundaryRegion copy in BoundaryWidth as a std::unique_ptr and I thought it would be nice to have unique_ptrs for the boundary and par_boundary vectors in Mesh/BoutMesh. I've implemented in 2e17480. I guess it breaks backward compatibility in principle, although I wouldn't expect end-users to need to call any of the changed functions, Mesh::getBoundaries() or Mesh::getBoundariesPar(). Should I keep or get rid of the commit?

• FieldData::setBoundary(region, op) was not implemented: in next at the moment it silently does nothing (except possibly printing "Replacing "). I've implemented what I think it was intended to do in ee22ccf but it might be better to just remove it. Obviously nobody is relying on it. Should I delete instead?

[johnomotani]

There was a comment in the manual about BoundaryWidth not being thread-safe because it changed the width of the BoundaryRegion while applying the BoundaryOp and then changed it back. That is no longer true, but are the boundary conditions thread-safe in general? I don't have a good understanding of what thread-safety means, but it looks to me like they are not, because they use shared BoundaryRegion objects as iterators, and BoundaryRegion::first(), BoundaryRegion::next1d(), etc. change the state of the shared object so if different threads were trying to apply boundary conditions to different fields at the same time, they might interfere with each other. Is that right? I've put a comment about the boundary conditions not being thread-safe in the manual, but someone who understands the concept better could probably write much more clearly... See https://github.com/boutproject/BOUT-dev/blob/ebda6f316a184e419876a573a4bdd58a44b425aa/manual/sphinx/user_docs/boundary_options.rst#limitations

[dschwoerer]

If you are not sure, could you not mention that at all?

Having wrong documentation is IMO way more confusing then no documentation. I don't think BOUT++ generally makes any sort of guarantees in that regard, thus saying one thing isn't thread safe implies the rest is - which might be the case, but I don't know.

[dschwoerer]

@johnomotani do you still intend to follow up on this?

[johnomotani]

@johnomotani do you still intend to follow up on this?

No. I don't remember now what exactly the purpose was, so don't know if in principle it would be worth updating and merging. In any case I don't have any time to work on it now, so I'm not going to do anything. Anyone is welcome to take it over if they think it's useful, although at this point it might be easier to start over from scratch!