Generalize the boundary condition options

[tberlok]

The boundary condition options are currently to fix variables to be zero at zmin and zmax or to not do anything.

This should be changed such that the boundary condition at zmin and zmax are specified independently. Furthermore, an option for specifying that the derivative of the variable is zero should be implemented.

[gkowal]

Thomas, I've started using psecas in my research and found it very useful. I need, however, to set both derivatives and values to zero at boundaries. Could you point me how to do this? Can I do this in solver._set_submatrix() by setting two first and two last columns to zero?

[tberlok]

Hi Grzegorz, great that you find psecas useful. Can you give me a few more details about what you need to do? That is, do you need to

1. set both the derivative and value of a variable to zero, or
2. either set the derivative or the value to zero.

If it is option 1, then I am not sure that I can quickly implement this functionality, if it is option 2 then I can probably try to implement it this week.

If it is option 1, and you are solving on a finite domain, then you might be able use the software dedalus instead of psecas (psecas was developed to handle problems on infinite and periodic domains). Alternatively, one could probably use this differentation matrix from dmsuite to modify psecas. I'd be interested in giving this a try but I don't know when I will have the time to do it.

[gkowal]

Thomas, in the problem I study I have an inflow at boundaries (Vz != 0 at z=-L and z=L). The domain is symmetric [-L,L]. This inflow varies depending on the solution of the interior. I believe at least the derivative of Vz should be set to zero, either on finite or infinite grid. Other variables should be zero at boundaries. The idea is to study how this inflow is affected by the solution near Z=0.

Thank you for pointing me to the Dedalus project. I will give it a try, but I order to help me understand better pseudospectral methods applied to my problem, I will try to incorporate the differentiation matrix from dmsuite in psecas.

[tberlok]

That sounds like an interesting problem. If only the derivative of Vz should be set to zero, then it sounds like option 2, which is what I intended to implement when I opened this Github issue. The method I had in mind is quite simple, but it requires modifying the right-hand side of the evp, turning it into a generalized evp. Those take longer to solve, so I like the idea of using cheb2c instead. Only caveat it that cheb2c only works on a finite domain.

The dmsuite differentation matrix is described in section 4 of Weidemann & Reddy where they also show how to solve a simple eigenvalue problem, their eq. 37. Awesome that you are interested in adding this to psecas. I took a quick look and wrote this down to help you get started:

from dmsuite import cheb2bc

# Example given by eq 37
# The boundary condition is u(-1) + u'(-1) = 0, u(1) = 0
# the point x=1 is simply removed, and xt therefore has N elements
N = 16
xt, d1t, d2t, phip, phim = cheb2bc(N, [[1, 1, 0], [1, 0, 0]])
e = np.linalg.eig(d2t)

# With boundary conditions u(-1) + u'(-1) = 0 and  u(1) + u'(1) = 0
# xt has N+1 elements, same as ChebyshevExtremaGrid
xt, d1t, d2t, phip, phim = cheb2bc(N, [[1, 1, 0], [1, 1, 0]])
from psecas import ChebyshevExtremaGrid
grid = ChebyshevExtremaGrid(N, zmin=-1, zmax=1)
# note that although xt and grid.zg contain the same values, xt is in reverse order.
# That is x[0]=1 while grid.zg[0] = -1

# Finally, for Neumann conditions, u'(-1) = 0 and  u'(1) = 0
# xt has N+1 elements, same as ChebyshevExtremaGrid

Depending on a given variable, a different differentation matrix would be used when evaluating its derivative. This would require making some changes to the equation parser in psecas. Let me know how it goes and if you have questions!

[tberlok]

Hi Grzegorz,

I have now made a quick implementation of my original idea. The code can now do Dirichlet or Neumann conditions. I have done minimal testing but the results look reasonable in my heat-flux-driven buoyancy instability example.

You can find the code in the branch neumann-boundary. The example is located at examples/galaxy-clusters/hbi.py and the underlying system class can be found at psecas/systems/hbi.py. Here I have set the boundary conditions using a new attribute, like this:

        # Boundary conditions
        self.boundaries = [False, True, True, False, True]

        # Extra information for boundary conditions
        self.extra_binfo = [[None, None], ['Dirichlet', 'Neumann'], 
                                     [None, 'Neumann'],
                            [None, None], ['Dirichlet', 'Neumann']]

This means that self.boundaries now just tells whether there is any boundaries to consider for a given variable. The details are then set in self.extra_info. The attribute self.boundaries still needs to be set and I have kept it for backwards compatibility (for now).

If you run the example examples/galaxy-clusters/hbi.py in a Python shell, then you'll find that

solver.system.result['dT'][0]  # is exactly zero
solver.system.result['dT'][-1] # is exactly zero
solver.system.result['dA'][0] # is exactly zero
grid.der(solver.system.result['dA'])[-1] # is very small 1e-12
solver.system.result['dvx'][0] # is not zero  (small though)
grid.der(solver.system.result['dvx'])[-1] # is very small 1e-12

which is what we would expect from the complicated boundary conditions.

This approach should also work on the infinite grids (I have not tested that yet, though) which I think is an advantage compared to the cheb2c approach.

One thing to note is that, for some problems, e.g. this one, setting too many or too few boundary conditions gives poor results. For instance, this example does not work well if a boundary condition is set on both dvx and dvz. So this is something to keep in mind.

I hope you can use this!

[gkowal]

Wow, Thomas, you are fast!

I didn't have chance to look into the code and think about what you suggested in your previous comment, and now you say you already implemented it. Great!

I will pull the branch neumann-boundary and test them with both finite and infinite grid.

Regarding your comment in the last paragraph. I've tested all combinations of True or False for each variable on infinite grid, and indeed some results simply didn't converge.

I'll give my feedback on Neumann boundaries soon.

[gkowal]

Thomas, I have tested my problem with both finite and infinite grid and Neumann boundary conditions. Although in the case of infinite grid the differences seem small, with the finite grid, I see clearly that the derivatives are zero but the values not. So, I think your quick implementation works well.

[tberlok]

That is great to hear! I hope you find some interesting results.

I also figured out how to add more non-trivial boundary conditions, e.g. yesterday I revisited an old problem on the domain r \in [1, 2] with the boundary condition r^2*dr(dr(u)) + r*dr(u) - u = 0 at r=1 and r=2. I will probably add an option to add such boundary strings. If you need this, let me know, and I will do it sooner rather than later.

If you have any other questions/problems you can always open an issue, or send me an email at my github username + gmail.com or aip.de.

[tberlok]

There is now an example of how to use boundary strings in /examples/accretion-disks/hall_mri.py.