Understanding the chebyshev solver in poisson2ND.py

[jchan192]

Hi,

I am learning the Chebyshev-Galerkin method in Shen95 for 2D, 3D Poisson Equation with nonhomogeneous Dirichlet BC. Your example poisson1D.py helps me a lot in understanding the method. But I am facing trouble in understanding poisson2ND.py, and hope that you can clarify something.

1. In the 2D case for Shen95, they are actually solving

In 1D, only the A matrix remains because alpha=0 for Poisson Eq. In both 2D and 3D cases, the B matrix is still there eventhough alpha=0. However, in your poisson2ND.py, you don't seem to have B matrix, but only used A by doing A = inner(v, -div(grad(u))) Is it true that we don't need B matrix, or it is already considered in your code somewhere that I do not understand?

2. Shen95 suggested using matrix decomposition, but la.solver2D seems to use LU factorization to solve the 2D poisson problem. In 1D Poisson problem can be easily solved using backward sweep method due to the tridiagonal property of A.
   Is it not possible to use backward sweeping method any more for 2D and 3D Poisson Eq.?

Thanks

[mikaem]

Hi
The assembly in poisson2ND returns a list of tensor product matrices and stores them in A. So A is just a name of a list, it is not a stiffness matrix in this demo. Each item in this list represents a term like BXB or AXB, where X is an outer (Kronecker) product and A and B are the matrices you refer to above. The list A also contains tensor product matrices for the boundary condition. The default solver la.Solver2D is a Kronecker product solver. That is, I assemble the complete 2D matrices from the 1D matrices A and B as Kronecker products. For this I'm using scipy.sparse.kron if I remember correctly. Hope this helps:-) You are more than welcome to implement a matrix decomposition solver. But since this is more work, and more difficult to generalize than the Kronecker product method, I have until now only implemented such a solver for a Legendre basis.

[jchan192]

Thanks for the reply and clarifying. I will certainly take some time to look at your matrix decomposition solver for the Legendre basis.

I have a question concerning about the BC in 2D. So in 1D, as mentioned an your document, the inhomogeneous BC can be easily assigned by u_hat[N-2] = a and u_hat[N-1] = b, after solving Au_hat = f_hat. I am wondering, conceptually, if the same procedure is also applicable to 2D case, so we just assign u_hat[N-2,:], u_hat[N-1,:], u_hat[:, N-2] and u_hat[:, N-1] by their corresponding BC, after solving Au_hatB + Bu_hatA^T = f_hat in 2D.

[mikaem]

@jchan192 It is more complicated than that in 2D. Conceptually you use a basis with two additional basis functions in your variational form for a given dimension. But then in 2D you use a basis that is the tensor product of 1D bases. This makes it more complicated. Unfortunately, I have not explained this very well anywhere. Closest thing is the Rayleigh-Benard demo, where inhomogeneous boundary conditions are coupled with a periodic Fourier space.