Nonuniform BoundaryConditions

[dschwoerer]

Add boundaries which are aware of non-uniform grid spacing:

• Dirichlet
• Neumann
• Free
• for 2nd, 3rd and 4th each

~Tests are not included, as they require realy~

The code uses sympy to generate the C++ code, so extend should be easy.

[johnomotani]

Definitely agree there is an issue here.

I think from a skim of your code here that the solution you implement is taking into account explicitly the positions of each point in x-y space. That, if I follow these lecture notes correctly --- https://www.nada.kth.se/kurser/kth/2D1263/l6.pdf [see section 4.3, in particular eqs (4.2), (4.4) and figs 3.3 and 3.5-6] --- is the most accurate thing to do, but our derivative operators don't do that, they just use the local grid spacing at the location of the derivative (updating the non-uniform grid support in the derivatives could be a useful thing to do...?).

Using the local spacing at the boundary is a bit of a pain, especially with staggered fields (they end up needing the grid spacing at the x-y corner of cells to put the x-boundary condition on a y-staggered field or vice versa). I have put together some changes, building on other stuff I was doing to get difops to converge accurately near boundaries, that get the grid spacings at the boundary locations. See the derivatives3 test in e8b28bbfed449077f34b723fbd47707b2809f8f3, which is updated to have a non-constant grid spacing and to apply boundary conditions. With the updates, it converges correctly for dirichlet, neumann and free boundary conditions, on staggered and unstaggered grids. The cost compared to your approach is having to change a lot more stuff in Coordinates, essentially to get a calcGridSpacing in src/mesh/boundary_op.cxx to be possible. (Although many of the changes were to enable the 'x-boundary of y-staggered field' case which should work but I haven't actually managed to test yet.) Some of the changes I think could be good anyway, from a neatness and/or performance viewpoint, but that's a question for a different PR (the commit history and tests definitely need tidying up).

[dschwoerer]

Some timings: - based on 41382ba6f21b13d725bf032a9de9af339aaea399

dirichlet_o2:   10000 iterations took 0.776276 s
neumann_o2:     10000 iterations took 2.3407 s
free_o2:        10000 iterations took 1.56801 s
dirichlet_nu_o2:        10000 iterations took 0.613323 s
neumann_nu_o2:  10000 iterations took 0.614806 s
free_nu_o2:     10000 iterations took 0.649469 s
dirichlet_o3:   10000 iterations took 5.79157 s
free_o3:        10000 iterations took 1.76174 s
dirichlet_nu_o3:        10000 iterations took 0.792008 s
neumann_nu_o3:  10000 iterations took 0.78357 s
free_nu_o3:     10000 iterations took 0.846108 s
dirichlet_o4:   10000 iterations took 6.96777 s
dirichlet_nu_o4:        10000 iterations took 1.02474 s
neumann_nu_o4:  10000 iterations took 1.54527 s
free_nu_o4:     10000 iterations took 1.06533 s

Only dirichlet_o4 is slightly faster then the more general versions. Thus ether the general versions could be improved (Is #1351 making them faster then the non-uniform versions?) or the more general version could be used as a default. If there is no performance penalty - I think we should try to do "the right thing"

I thought of implementing some adaptive stencils - but haven't came around to do it. I thought of doing it as an aiolos-only feature, as the standard stencil code does not know the position of the derivative.

[dschwoerer]

This is ready to be merged - can someone have a look over this?

As this is to my understanding the correct implementation, thus we could discuss to replace the old, non-uniform mesh aware implementation, but I thought as a first step we could have it in parallel, to allow testing ...

[ZedThree]

It would be really good if this included the tests as well. Can they be written without using realy?

[johnomotani]

@dschwoerer could you add a docstring explaining how the coefficients are calculated? I'd find that very helpful.

[dschwoerer]

Do you mean the derivation, starting with the taylor expansion, or rather explain the implementation?

[johnomotani]

If you could sketch the derivation, with notation corresponding to the variables in the code, that would be great :1st_place_medal:

[dschwoerer]

I think I should have addressed all comments. Sorry for the delay, but I stopped using non-uniform meshes for a while ...

[dschwoerer]

Can someone give the failed build a kick?

@johnomotani is the description in bfc4cbc sufficient?

[dschwoerer]

Boundary condition	1000 ; 2D ; 128^3	1000 ; 2D ; 128x16384x1	1000 ; 3D ; 128^3
dirichlet_o2	0.353929 s	0.476291 s	0.375638 s
neumann_o2	0.442191 s	0.5878 s	0.856184 s
free_o2	0.217463 s	0.421346 s	0.243244 s
dirichlet_nu_o2	0.187466 s	1.97744 s	2.78296 s
neumann_nu_o2	0.177839 s	1.79646 s	2.10596 s
free_nu_o2	0.187784 s	2.74089 s	3.44364 s
dirichlet_o3	0.829212 s	0.785762 s	0.907647 s
free_o3	0.309709 s	0.984481 s	0.382664 s
dirichlet_nu_o3	0.251681 s	3.85593 s	5.58085 s
neumann_nu_o3	0.238788 s	2.94091 s	3.49203 s
free_nu_o3	0.287577 s	3.22091 s	4.60549 s
dirichlet_o4	1.09995 s	0.773845 s	1.07976 s
dirichlet_nu_o4	0.339356 s	3.487 s	4.87303 s
neumann_nu_o4	0.362675 s	4.04734 s	6.26752 s
free_nu_o4	0.386731 s	3.74004 s	5.49017 s

So for 2D and large enough nz the overhead of the correct stencil is minimal. Otherwise it can be expensive. Would be nice to know why 3D is slower, probably something that can be fixed.