spectralDNS / shenfun

High performance computational platform in Python for the spectral Galerkin method
http://shenfun.readthedocs.org
BSD 2-Clause "Simplified" License
204 stars 42 forks source link

Old shenfun script for 2D Poisson equation dosn't work with new version #83

Closed yesint closed 3 years ago

yesint commented 3 years ago

I apologize for possible stupid question. I'm not familiar to the theory under shenfun at all and using it as a black box. Previously I was using the following script to solve 2D Poisson equation with periodic boundary conditions (data preparation part omitted):

K1 = FunctionSpace(N[0], family='F', dtype='D',  domain=(0, Xsz))
K2 = FunctionSpace(N[1], family='F', dtype='d', domain=(0, Zsz))
T = TensorProductSpace(comm, (K1, K2), axes=(0,1))
X = T.local_mesh(True)
u = TrialFunction(T)
v = TestFunction(T)

# Get f
fj = Array(T, buffer=fe)

# Compute right hand side of Poisson equation
f_hat = Function(T)
f_hat = inner(v, fj, output_array=f_hat)

# Solve Poisson equation
u_hat = Function(T)
A = inner(grad(v), grad(u))
u_hat = A.solve(-f_hat, u_hat)

uq = Array(T)
uq = T.backward(u_hat, uq, fast_transform=True)

it was working perfectly fine before, but with new version of Shenfun it fails:

Traceback (most recent call last):
  File "shen2d-2.py", line 55, in <module>
    u_hat = A.solve(-f_hat, u_hat)
AttributeError: 'list' object has no attribute 'solve'

How can I fix it?

mikaem commented 3 years ago

Yes, the linear algebra solvers have changed interface lately. I'm trying to keep the demos up-to-date, so you should have a look at the 2D Poisson demo for an update. Basically, use

A = inner(grad(v), grad(u))
sol = la.SolverDiagonal(A)
u_hat = sol(-f_hat, u_hat, constraints=((0, 0),))

Alternatively you can do

from shenfun import get_simplified_tpmatrices
B = get_simplified_tpmatrices(A)[0]
u_hat = B.solve(-f_hat, u_hat)

Before we used to return simplified matrices from inner. Simplified extracts the diagonals of the Fourier matrices and places them in scale arrays in order to facilitate a fast inverse solve, which is then basically just f_hat/B.scale, taking care of wavenumber 0. But I guess you need to understand how shenfun uses 1D matrices and tensor product matrices (TPMatrix) in order to fully make sense of this:-)

yesint commented 3 years ago

Thanks a lot! It works.