spectralDNS / shenfun

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

modifying dirichlet_poisson2D.py #42

Open francispoulin opened 4 years ago

francispoulin commented 4 years ago

I am modifying dirichlet_poisson2D.py to solve a different example and I'm having some difficulties. Maybe it has to do with the fact that my right-hand-side, f, is independent of one of the two variables.

Below is the error message that I obtain as well as the parts fo the code that I modified, with the lines that I have removed simply commented out.

I suspect I am making a rookie mistake here but cannot yet figure out what that is. Could you help me understand what I am doing wrong here?


`#Error message

$ python dirichlet_poisson2D.py 32 chebyshev
Traceback (most recent call last):
  File "dirichlet_poisson2D.py", line 58, in <module>
    fj = Array(T, buffer=fe)
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/shenfun/forms/arguments.py", line 816, in __new__
    buffer = sympy.lambdify(sym0, buffer)(*space.local_mesh())
TypeError: _lambdifygenerated() takes 1 positional argument but 2 were given

`

`# snipit of code

a = -1.
b =  1.
x, y = symbols("x,y")
ue = tanh(x)  
#ue =  (cos(2.*np.pi*x/10.) + sin(y))*(1 - x**2) - a*(1 - x)/2. - b*(1 + x)/2.
#fe = ue.diff(x, 2) + ue.diff(y, 2)
fe = -2.*tanh(x)/pow(cosh(x),2)

# Size of discretization
N = (int(sys.argv[-2]), int(sys.argv[-2])+1)

#SD = Basis(N[0], family=family, scaled=True, bc=(a, b))
SD = Basis(N[0], family=family, domain=(-10., 10.), scaled=True, bc=(a, b))
K1 = Basis(N[1], family='F', dtype='d', domain=(0., 2*np.pi))
T = TensorProductSpace(comm, (SD, K1), axes=(0, 1))
X = T.local_mesh(True)
u = TrialFunction(T)
v = TestFunction(T)

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

`

mikaem commented 4 years ago

Yes, you're absolutely right, but I', pretty sure this has already been fixed, so are you on an old version of shenfun?

francispoulin commented 4 years ago

Ah, I wasn't aware of the update.

I just updates shenfun using conda and I still obtain the same error.

Below is a list of the packages that I believe were updated. Is there one missing that should fix this problem?

`# Updated packages

The following packages will be UPDATED:

  certifi                                 2019.11.28-py38_0 --> 2020.4.5.1-py38_0
  curl                                    7.68.0-hbc83047_0 --> 7.69.1-hbc83047_0
  libcurl                                 7.68.0-h20c2e04_0 --> 7.69.1-h20c2e04_0
  libssh2                                  1.8.2-h1ba5d50_0 --> 1.9.0-h1ba5d50_1
  openssl                                 1.1.1d-h7b6447c_4 --> 1.1.1g-h7b6447c_0
  packaging                                       20.1-py_0 --> 20.3-py_0
  python                                   3.8.1-h0371630_1 --> 3.8.2-hcf32534_0
  readline                                   7.0-h7b6447c_5 --> 8.0-h7b6447c_0
  setuptools                                  45.2.0-py38_0 --> 46.1.3-py38_0
  sqlite                                  3.31.1-h7b6447c_0 --> 3.31.1-h62c20be_1
  tornado                              6.0.3-py38h7b6447c_3 --> 6.0.4-py38h7b6447c_1
  wcwidth                                        0.1.8-py_0 --> 0.1.9-py_0
  xz                                       5.2.4-h14c3975_4 --> 5.2.5-h7b6447c_0

The following packages will be SUPERSEDED by a higher-priority channel:

kiwisolver conda-forge::kiwisolver-1.1.0-py38hc9~ --> pkgs/main::kiwisolver-1.0.1-py38he6710b0_0 `

francispoulin commented 4 years ago

Incidently, I did get this working in the 1D version so happy to learn how to do that correctly.

mikaem commented 4 years ago

Ok, then I need to update the conda package:-) It is fixed in master. Meanwhile, just use something with 0y (or 1e-16y) to get it to work:-)

francispoulin commented 4 years ago

Thanks. I did try that, see below, and that did not fix the error unfortunately.

`# My attempt

ue = tanh(x) + 0*y

`

mikaem commented 4 years ago

You did not try 1e-16y😀

mikaem commented 4 years ago

Sympy is smart enough to eliminate the term that is identically 0.

francispoulin commented 4 years ago

I just tried now, sorry for not trying it before, and the problem persisted. Even when I increased it to 1e-8*y the same error occurred.

mikaem commented 4 years ago

ue is differentiated twice😃

francispoulin commented 4 years ago

Right, of course, so that term will disappear.

In this version, I didn't compute the derivative but specified the fe exactly.

I am happy to say that when I add 1e-16 to both ue and fe it does work. I will look at the convergence to make sure that everything looks good.

I'll keep this open for the moment but will close it when I get the update.

Thanks again!