2D Rayleigh-Benard in vorticity-streamfunction formulation

[loiseaujc]

Hi Mikael.

Based on shenfun, I have been able to implement a fairly efficient solver for two-dimensional decaying turbulence in a doubly-periodic domain which I eventually plan to use for part of my lectures. I am now willing to use shenfun again to implement a two-dimensional Rayleigh-Bénard solver using a vorticity-streamfunction formulation and I have a few questions. Before getting there, here are the equations I am interested in :

where is the out-of-plane vorticity, is the streamfunction and is the deviation from the linear temperature profile. The set of boundary conditions needed to close the system are the following :

• Homogeneous dirichlet conditions for on both walls
• Bihamonic boundary conditions for on both walls as well.

Finally, all variables are periodic in the horizontal direction. The problem we consider is thus the canonical two-dimensional Rayleigh-Bénard convection with solid isothermal walls.

As far as my understanding of shenfun goes, here is basically what I need to do then :

• I can use the same Fourier basis for all three variables, i.e.
  F = Basis(N[1], family='F', dtype='D', domain=(-np.pi, np.pi))

• I then need to define the appropriate Chebyshev basis for each variable as
  C_vorticity = Basis(N[0], family='C', domain=(0, 1))
C_temperature = Basis(N[0], family='C', domain=(0, 1), bc='dirichlet')
C_stream = Basis(N[0], family='C', domain=(0, 1), bc='biharmonic')

• Given these different bases, I can then define the appropriate TensorProductSpaces as
  V_vorticity = TensorProductSpace(comm, (C_vorticity, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))
V_temperature = TensorProductSpace(comm, (C_temperature, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))
V_stream = TensorProductSpace(comm, (C_stream, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))

• Finally, I can create the VectorProductSpace for my state vector as
  Q = VectorTensorProductSpace([V_vorticity, V_temperature])

If all of this is correct, I then have two questions, one related to the LinearRHS function, and the other to the NonlinearRHS ones. First, about LinearRHS : given my VectorTensorProductSpace Q, how can I create the linear operator

( where I assume that I treat and explicitly in the NonlinearRHS function)?

If creating such a block-diagonal linear operator is currently doable in shenfun, which type of linear solver does shenfun then uses to invert it considering the two different TensorProductSpaces used?

Thanks a lot,
JC

EDIT : I may have mixed some signs, but you get the points anyway.

[mikaem]

Dear JC,

Very cool that you are implementing this using shenfun:-) Just so as you know, I recently added two very efficient 3D Rayleigh Benard solvers here and here. A demo running these solvers is here

Regarding your questions about LinearRHS and NonlinearRHS. Unfortunately, the few implemented integrators in Shenfun are only valid for pure Fourier problems, since for Dirichlet comes the added complexity of a non-diagonal mass matrix. With your equations you basically have Helmholtz and biharmonic equations to be solved for each time step and there are very efficient (O(N)) solvers available for this in shenfun. But you need to set up the equations yourselves, e.g., using Crank-Nicolson, backwards differencing or similar. Not very difficult, I have done it for the channel solver (KMM) and the 3D Rayleigh Benard solvers, but there is some book-keeping and everything is not as simple as for the steady demo programs:-) Solving two equations coupled is also something I have been wanting to do, but currently there is little there to help, sorry:-(

[loiseaujc]

I was expecting something along these lines (particularly about the integrators), but thanks anyway for the rapid answer.

I will work out the weak formulation for the linearized equations first tomorrow as I get to the lab, and implement the linearized solver just to make sure everything works out and I will see for the nonlinear term later on.

As for the solver you have already implemented, I already took a look at them and they seem pretty cool indeed. The reason I want to write my own version is to make it a bit more readable for students (or colleagues) not as well versed into python as we can be :-]

[mikaem]

Great:-) A 2D R-B with vorticity-streamfunction formulation would make for a great demo I think. A time-dependent inhomogeneous flow is currently sorely missed. The spectralDNS solvers are in deed to complex, because all bits and pieces are torn apart for speed.

[loiseaujc]

Hi Mikael,

I have started to work on the linearized Rayleigh-Benard equations. Using Crank-Nicholson for the diffusive terms and Adam-Bashforth for the coupling ones, here is my weak formulation for the streamfunction equation:

or equivalently in compact form

I see how to construct each matrix independently:

• and are simply two instances of the Biharmonic class where I change the scale attribute attached to each submatrices accordingly.
• is constructed as inner( v_psi , Derivative(u_theta, 0, k=1) ).

Inverting and computing the matrix-vector product involving is quite straightforward. My problem is with the computation of the matrix-vector product involving . I dug into the sources of shenfun and it seems that Biharmonic has no matvec method. For the evolution equation of the temperature, while Helhmholtz does have a matvec method, it seems to be implemented only for one-dimensional fields, but not for 2D or 3D ones. Is there any reason for that or did I miss something?

I could always use the lower-level functions but I was willing to try to keep the code as high-level as possible for readability.

[mikaem]

Hi
You are correct. The matrix-vector product of the biharmonic system has not made it all into the repository yet. But you are in luck, because I'm currently working on moving the channel flow solver in spectralDNS such that the inhomogeneous direction is in the last, and not first, axis. There is a 1D and 3D biharmonic matvec version in spectralDNS. I plan to implement this in shenfun and to make it generic in terms of axes and dimension. Should be done today. Same for Helmholtz I guess.

However, note that all matrices in the system, the biharmonic, stiffness and mass matrix, they all should have fast implementations of the matvec. So you could split it up and do each term by itself. See SBBmat, ABBmat and BBBmat in chebyshev.matrices.py. Hmm. I now see that fast implementation is lacking for 2D. So that needs to be done as well:-)

[mikaem]

I've pushed matvecs to both Biharmonic and Helmholtz solvers. They're not all optimal, but it should work for all axes and dimensions.

[loiseaujc]

Awesome! I will actually fork and install shenfun locally rather than going through conda.
I was starting to implement my own work-around in the mean time but you've been much faster than I did :] I will give it a try right away.

[mikaem]

Just a quick fix in Legendre. But these matrices are very sparse and all have highly efficient matvecs anyway. Just the Chebyshev matrices that require some effort really.

[loiseaujc]

Either I am missing something, or there might be a small issue with the matrix-vector product. Starting from the 2D Poisson-Dirichlet example, here is what I did to test matvec:

In [11]: f_bis = f_hat.copy()

In [12]: f_bis.fill(0)

In [13]: f_bis = H.matvec(u_hat, f_bis)

In [14]: abs(f_bis-f_hat).max()
Out[14]: Function(7.85398163)

I would expect f_bis and f_hat to be equal. When I plot the difference between these two arrays in physical space, it looks like delta(x, y) = x (where x is the chebyshev direction). Note that I have the same results if I compute the corresponding matrix-vector product directly from the lower-level methods attached to each type of matrices.

EDIT: In Ipython, my very first command is obviously run dirichlet_poisson2D.py 32 chebyshev. Note that I still have the problem if I use Legendre polynomials instead.

[mikaem]

Ok. It's probably due to inhomogeneous boundary conditions. Those require adjustments that I did not consider yet. Could you please test with zero bcs? Don't have my computer.

[loiseaujc]

I ran the same test with the 1D Poisson solver and homogeneous dirichlet boundary conditions. Here is the output:

In [72]: run dirichlet_poisson1D.py 16 chebyshev
Error=5.2901900972253749e-06

In [73]: f_bis = f_hat.copy()

In [74]: f_bis.fill(0)

In [75]: f_bis = A.matvec(u_hat, f_bis)

In [76]: f_bis-f_hat
Out[76]: 
Function([ 0.00000000e+00,  0.00000000e+00, -1.57772181e-30,
           0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00, -1.57772181e-30,
           0.00000000e+00, -9.86076132e-32,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00, -1.31703830e+01,
           1.31703830e+01])

It indeed seems like the problem comes from the boundary conditions, whether they are homogeneous or not.

[mikaem]

Yes. The last two items are from bcs. There should be no error for a homogeneous case. Did you try that?

[loiseaujc]

That test was for homogeneous boundary conditions. I modified the example as

SD = Basis(N, family=family, domain=domain, bc='dirichlet')

Basically, the last two entries of f_bis = A.matvec(u_hat, f_bis) are set to 0 (which I guess is consistent with the boundary conditions), while the Chebyshev transform f_hat is simply the Chebyshev transform. Does that make sense? If that is the case, then it ain't a problem since, whenever I will solve the Helmholtz equation in my time-stepping, the last two entries of the right-hand side vector will be set to whatever is needed for the boundary conditions for the updated velocity. Is that correct?

[mikaem]

The thing is, it's not really an error. Those last two items of the f_hat array are not in use when solving a problem with inhomogeneous bcs. Because those two items are set equal to a and b during solve. They are not considered unknowns. I haven't thought about this regarding matvecs yet, but it should be easy enough to fix. Just disregard the last two items in the Dirichlet direction for now.