spectraldns / shenfun Goto Github PK

View Code? Open in Web Editor NEW

184.0 14.0 42.0 93 MB

High performance computational platform in Python for the spectral Galerkin method

Home Page: http://shenfun.readthedocs.org

License: BSD 2-Clause "Simplified" License

Python 66.72% Shell 0.06% Makefile 0.02% Jupyter Notebook 11.78% C++ 9.61% Cython 11.81%

spectral-methods galerkin fenics turbulence fourier spectral python chebyshev legendre laguerre

shenfun's Introduction

Shenfun

Try it in a jupyter hub using Binder

Description

Shenfun is a high performance computing platform for solving partial differential equations (PDEs) by the spectral Galerkin method. The user interface to shenfun is very similar to FEniCS, but applications are limited to multidimensional tensor product grids, using either Cartesian or curvilinear grids (e.g., but not limited to, polar, cylindrical, spherical or parabolic). The code is parallelized with MPI through the mpi4py-fft package.

Shenfun enables fast development of efficient and accurate PDE solvers (spectral order and accuracy), in the comfortable high-level Python language. The spectral accuracy is ensured by using high-order global orthogonal basis functions (Fourier, Legendre, Chebyshev first and second kind, Ultraspherical, Jacobi, Laguerre and Hermite), as opposed to finite element codes that are using low-order local basis functions. Efficiency is ensured through vectorization (Numpy), parallelization (mpi4py) and by moving critical routines to Cython or Numba. Shenfun has been used to run turbulence simulations (Direct Numerical Simulations) on thousands of processors on high-performance supercomputers, see the spectralDNS repository.

The demo folder contains several examples for the Poisson, Helmholtz and Biharmonic equations. For extended documentation and installation instructions see ReadTheDocs. For interactive demos, see the jupyter book. Note that shenfun currently comes with the possibility to use two non-periodic directions (see biharmonic demo), and equations may be solved coupled and implicit (see MixedPoisson.py).

Note that shenfun works with curvilinear coordinates. For example, it is possible to solve equations on a sphere (using spherical coordinates), on the surface of a torus, on a Möbius strip or along any curved line in 2D/3D. Actually, any new coordinates may be defined by the user as long as the coordinates lead to a system of equations with separable coefficients. After defining new coordinates, operators like div, grad and curl work automatically with the new curvilinear coordinates. See also this notebook on the sphere or an illustration of the vector Laplacian.

For a more psychedelic experience, have a look at the simulation of the Ginzburg-Landau equation on the sphere (click for Youtube-video):

Shenfun can also be used to approximate analytical functions with global spectral basis functions, and to integrate over highly complex domains, like the seashell below, see this demo.

Installation

Shenfun can be installed using either pip or conda, see installation chapter on readthedocs.

Dependencies

Python 3.7 or above. Test suits are run with Python 3.8, 3.9 and 3.10.

A functional MPI 2.x/3.x implementation like MPICH or Open MPI built with shared/dynamic libraries.

FFTW version 3, also built with shared/dynamic libraries.

Python modules:

Numpy

Scipy

Sympy

Cython

mpi4py

mpi4py-fft

Contact

For comments, issues, bug-reports and requests, please use the issue tracker of the current repository, or see How to contribute? at readthedocs. Otherwise the principal author can be reached at:

Mikael Mortensen
mikaem at math.uio.no
https://mikaem.github.io/
Department of Mathematics
University of Oslo
Norway

shenfun's People

Contributors

Stargazers

Watchers

shenfun's Issues

Solver for linear system of a VectorTensorProductSpace

Dear Mikael,

I would like to solve a vector Poisson problem in 2D. For now, I am solving the vector Poisson problem as a test case, where the spatial components do not couple. Later, I would like to move on to other problems, in which the spatial components are coupled.

I am using Legendre polynomials and I have created the standard weak form using integration by parts: inner(grad(u), grad(v)) where the functions u and v are elements of a VectorTensorProductSpace. As far as I understand, the function inner is used to compute the components of the linear system. In some of the demos, specialized classes like BlockMatrix are used to create the linear system which is finally solved. However, the BlockMatrix class only seems to work for a MixedTensorProductSpace. In my case, the inner function returns a list of TPMatrix objects and I am not sure what exactly they represent. Moreover, it is not clear how to compose a linear system from them. I hope my code below shows the problem.

Best regards,
Thilo

from mpi4py import MPI
from sympy import symbols
from shenfun import inner, grad, TestFunction, TrialFunction, Array, \
    Basis, TensorProductSpace, VectorTensorProductSpace
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD
family = 'legendre'

# define rhs
x, y = symbols("x,y")
f = (x**2, y)

# Size of discretization
n = 25
N = [n, n+1]

# Create bases and spaces
BX = Basis(N[0], family=family, bc=(0, 0))
BY = Basis(N[1], family=family, bc=(0, 0))

subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T = TensorProductSpace(subcomms, (BX, BY), axes=(1, 0))
V = VectorTensorProductSpace([T, T])

u = TrialFunction(V)
v = TestFunction(V)

# get f on quad points
fj = Array(V, buffer=f)

# compute matrices and rhs
matrices = inner(grad(v), grad(u))
b = inner(v, fj)

spectralDNS tests fail because it expects the function pyfftw.builders.dct

When I run pytest -v I see a lot of test failures (see test summary below):

test_MHD.py::test_MHD[uniform] PASSED                                                                                                                                                          [  2%]
test_MHD.py::test_MHD[nonuniform] PASSED                                                                                                                                                       [  4%]
test_NS2D.py::test_NS2D[1] PASSED                                                                                                                                                              [  6%]
test_NS2D.py::test_NS2D[2] PASSED                                                                                                                                                              [  8%]
test_NSVV.py::test_solvers[NS_uniform] PASSED                                                                                                                                                  [ 10%]
test_NSVV.py::test_solvers[VV_uniform] PASSED                                                                                                                                                  [ 12%]
test_NSVV.py::test_solvers[NS_nonuniform] PASSED                                                                                                                                               [ 14%]
test_NSVV.py::test_solvers[VV_nonuniform] PASSED                                                                                                                                               [ 16%]
test_NSVV.py::test_integrators[NS_uniform] FAILED                                                                                                                                              [ 18%]
test_NSVV.py::test_integrators[VV_uniform] FAILED                                                                                                                                              [ 20%]
test_NSVV.py::test_integrators[NS_nonuniform] FAILED                                                                                                                                           [ 22%]
test_NSVV.py::test_integrators[VV_nonuniform] FAILED                                                                                                                                           [ 24%]
test_channel.py::test_channel[KMMRK3] FAILED                                                                                                                                                   [ 26%]
test_channel.py::test_channel[KMM] FAILED                                                                                                                                                      [ 28%]
test_channel.py::test_channel[IPCS] SKIPPED                                                                                                                                                    [ 30%]
test_channel.py::test_channel[IPCSR] SKIPPED                                                                                                                                                   [ 32%]
test_channel.py::test_channel[KMMRK3_mpifft4py] FAILED                                                                                                                                         [ 34%]
test_channel.py::test_channel[KMM_mpifft4py] FAILED                                                                                                                                            [ 36%]
test_shentransforms.py::test_FST[GC-Basis] FAILED                                                                                                                                              [ 38%]
test_shentransforms.py::test_FST[GC-ShenDirichletBasis] FAILED                                                                                                                                 [ 40%]
test_shentransforms.py::test_FST[GC-ShenNeumannBasis] FAILED                                                                                                                                   [ 42%]
test_shentransforms.py::test_FST[GC-ShenBiharmonicBasis] FAILED                                                                                                                                [ 44%]
test_shentransforms.py::test_FST[GL-Basis] FAILED                                                                                                                                              [ 46%]
test_shentransforms.py::test_FST[GL-ShenDirichletBasis] FAILED                                                                                                                                 [ 48%]
test_shentransforms.py::test_FST[GL-ShenNeumannBasis] FAILED                                                                                                                                   [ 51%]
test_shentransforms.py::test_FST[GL-ShenBiharmonicBasis] FAILED                                                                                                                                [ 53%]
test_shentransforms.py::test_FST_padded[GC-Basis] FAILED                                                                                                                                       [ 55%]
test_shentransforms.py::test_FST_padded[GC-ShenDirichletBasis] FAILED                                                                                                                          [ 57%]
test_shentransforms.py::test_FST_padded[GC-ShenNeumannBasis] FAILED                                                                                                                            [ 59%]
test_shentransforms.py::test_FST_padded[GC-ShenBiharmonicBasis] FAILED                                                                                                                         [ 61%]
test_shentransforms.py::test_FST_padded[GL-Basis] FAILED                                                                                                                                       [ 63%]
test_shentransforms.py::test_FST_padded[GL-ShenDirichletBasis] FAILED                                                                                                                          [ 65%]
test_shentransforms.py::test_FST_padded[GL-ShenNeumannBasis] FAILED                                                                                                                            [ 67%]
test_shentransforms.py::test_FST_padded[GL-ShenBiharmonicBasis] FAILED                                                                                                                         [ 69%]
test_shentransforms.py::test_Mult_Div FAILED                                                                                                                                                   [ 71%]
test_shentransforms.py::test_Helmholtz[GC-ShenDirichletBasis] FAILED                                                                                                                           [ 73%]
test_shentransforms.py::test_Helmholtz[GC-ShenNeumannBasis] FAILED                                                                                                                             [ 75%]
test_shentransforms.py::test_Helmholtz[GL-ShenDirichletBasis] FAILED                                                                                                                           [ 77%]
test_shentransforms.py::test_Helmholtz[GL-ShenNeumannBasis] FAILED                                                                                                                             [ 79%]
test_shentransforms.py::test_Helmholtz2[GC] FAILED                                                                                                                                             [ 81%]
test_shentransforms.py::test_Helmholtz2[GL] FAILED                                                                                                                                             [ 83%]
test_shentransforms.py::test_Mult_CTD[GC] FAILED                                                                                                                                               [ 85%]
test_shentransforms.py::test_Mult_CTD[GL] FAILED                                                                                                                                               [ 87%]
test_shentransforms.py::test_Mult_CTD_3D[GC] FAILED                                                                                                                                            [ 89%]
test_shentransforms.py::test_Mult_CTD_3D[GL] FAILED                                                                                                                                            [ 91%]
test_shentransforms.py::test_Biharmonic[GC] FAILED                                                                                                                                             [ 93%]
test_shentransforms.py::test_Biharmonic[GL] FAILED                                                                                                                                             [ 95%]
test_shentransforms.py::test_Helmholtz_matvec[GC] FAILED                                                                                                                                       [ 97%]
test_shentransforms.py::test_Helmholtz_matvec[GL] FAILED                                                                                                                                       [100%]

Most of the failures are due to AttributeErrors like:

self = <shenfun.chebyshev.bases.Basis object at 0x7fd643bd4910>, N = 64, quad = 'GL', plan = True, domain = (-1.0, 1.0)

    def __init__(self, N=0, quad="GC", plan=False, domain=(-1., 1.)):
        ChebyshevBase.__init__(self, N, quad, domain)
        if quad == 'GC':
            self._xfftn_fwd = functools.partial(pyfftw.builders.dct, type=2)
            self._xfftn_bck = functools.partial(pyfftw.builders.dct, type=3)
        else:
>           self._xfftn_fwd = functools.partial(pyfftw.builders.dct, type=1)
E           AttributeError: 'module' object has no attribute 'dct'

/scratch/avmo/opt/spectraldns/lib/python2.7/site-packages/shenfun/chebyshev/bases.py:202: AttributeError

I am using the latest version of pyfftw

» pip show pyfftw
Name: pyFFTW
Version: 0.10.4
Summary: A pythonic wrapper around FFTW, the FFT library, presenting a unified interface for all the supported transforms.
Home-page: http://hgomersall.github.com/pyFFTW/
Author: Henry Gomersall
Author-email: [email protected]
License: UNKNOWN
Location: /scratch/avmo/opt/spectraldns/lib/python2.7/site-packages
Requires:

And even in their github master version I don't see this dct function: https://github.com/pyFFTW/pyFFTW/blob/master/pyfftw/builders/builders.py. Are you using a modified version of pyFFTW?

Two inhomogeneous directions?

Hej Mikael,

I have successfully written as small python code for two-dimensional turbulence using shenfun. I am planning to move on to a fully three-dimensional channel flow, and maybe eventually a duct flow. In that matter, it is written in the abstract of the shenfun paper:

With the shenfun Python module (github.com/spectralDNS/shenfun) an effort is made towards automating the implementation of the spectral Galerkin method for simple tensor product domains, consisting of (currently) one non-periodic and any number of periodic directions.

Do you have any plan of adding eventually multiple (say just two to begin with) inhomogeneous directions? That would be extremely valuable for people interested in confined flows (e.g. lid-driven cavity flows or duct flow to mention just two of them).

Thanks a lot anyway for developing shenfun. A really cool package!

Problem with Chebyshev basis?

Hi Mikael,

I have started to implement my Rayleigh-Benard solver and wanted to play a bit with the Poisson and Biharmonic solvers first. To do so, I downloaded the scripts available in the demo folder and I have some weird things going on:

For both the Poisson and Biharmonic scipts (2D ones), if I set the basis used in the non-homogeneous direction to be Chebyshev, the right-hand-side vector f_hat = inner(v, fj, output_array=f_hat) all its entries turn out be to equal to 0 + 0j.
If instead I set the basis to be Legendre polynomials, then everything works out correctly.

For information, I have installed shenfun using conda.

improving scalability

Hello Mikael,

I have tested the scalability of my 2D Navier-Stokes (shared in another issue) code and find that it's not doing great. It seems to get worst as the number of degrees of freedom increases. I have copied the results in a table below in an text table.

Just to be clear, I wasn't experting great efficiency as I have not been very clever in the setup of the code, but I do find it strange that as the number of points increase, the efficiency decrease. Normally, I would expect it to get better.

Question 1: Does this look odd to you?

Question 2: If I wanted to monitor the efficiency of my code in parallel, what would you recommend I use?

Note: these are done on my new, person desktop with 18 cores and 64 GB of ram. I am trying to get shenfun installed on a server, which has more cores to play with, but the people in charge have decided not to support conda, and don't want us to install conda, so that is becoming more of a chore than it should be. This may make it's way into another issue sometime ;)

| Np=1 | 2 | 4 | 8 | 16 |
N=1024 | 1 | 0.68 | 0.61 | 0.4 | 0.26
N=2048 | 1 | 0.55 | 0.47 | 0.37 | 0.19
N=4096 | 1 | 0.47 | 0.39 | 0.27 | 0.16

different results in parallel + error in output

Maybe I should have split this up but both happen when I go from one course to two, so maybe they are connected.

First, the results in serial look great and are way faster than was I was doing previously. I think it goes to show that RK4 is much faster than Adams Bashforth 3.

Unfortunately, when I run it in parallel i get two odd results.

First, the writing seems to freeze and I need to kill it. The results are copied below.

Second, I am looking at a double jet and their initial conditions are almost identical except for a random perturbation, but as one can see from a snapshot they evolve differently. I would say one is more damped then the other, which might suggest something to do with the padding. But I tried it without padding and it did not resolve the issue. Unfortunately, it's not clear to me how to include a figure in this post but I am happy to do that would help.

Clearly I've messed something up but not sure what. Any suggestions what I should look into?

`# error

/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/mpi4py_fft/io/generate_xdmf.py:117: H5pyDeprecationWarning: The default file mode will change to 'r' (read-only) in h5py 3.0. To suppress this warning, pass the mode you need to h5py.File(), or set the global default h5.get_config().default_file_mode, or set the environment variable H5PY_DEFAULT_READONLY=1. Available modes are: 'r', 'r+', 'w', 'w-'/'x', 'a'. See the docs for details.
f = h5py.File(h5filename)
Traceback (most recent call last):
File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/h5py/_hl/files.py", line 199, in make_fid
fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 88, in h5py.h5f.open
OSError: Unable to open file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/h5py/_hl/files.py", line 202, in make_fid
fid = h5f.open(name, h5f.ACC_RDONLY, fapl=fapl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 88, in h5py.h5f.open
OSError: Unable to open file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "NavierStokes2D.py", line 141, in
fftw.export_wisdom('QG1L.wisdom')
File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/mpi4py_fft/io/generate_xdmf.py", line 117, in generate_xdmf
f = h5py.File(h5filename)
File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/h5py/_hl/files.py", line 406, in init
fid = make_fid(name, mode, userblock_size,
File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/h5py/_hl/files.py", line 204, in make_fid
fid = h5f.create(name, h5f.ACC_EXCL, fapl=fapl, fcpl=fcpl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to open file: name = 'QG_1Layer_256.h5', errno = 17, error message = 'File exists', flags = 15, o_flags = c2)
`

10 errors in tests

I have installed shenfun using conda, as requested and that seemed to go well. However, when I try the tests, as suggested I get 10 errors. Maybe my conda install wasn't as good as I thought it was?

`$ python -m pytest tests/
=============================== test session starts ===============================
platform linux -- Python 3.8.2, pytest-5.3.5, py-1.8.1, pluggy-0.13.1
rootdir: /home/fpoulin/software/shenfun
collected 1 item / 10 errors

===================================== ERRORS ======================================
______________________ ERROR collecting tests/test_bases.py _______________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_bases.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_bases.py:3: in
from shenfun import *
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
_______________________ ERROR collecting tests/test_curl.py _______________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_curl.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_curl.py:4: in
from shenfun import inner, div, curl, TestFunction, TrialFunction, Function,
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
______________________ ERROR collecting tests/test_forms.py _______________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_forms.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_forms.py:4: in
import shenfun
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
_____________________ ERROR collecting tests/test_fourier.py ______________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_fourier.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_fourier.py:2: in
from shenfun import *
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
________________________ ERROR collecting tests/test_io.py ________________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_io.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_io.py:6: in
from shenfun import *
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
________________________ ERROR collecting tests/test_la.py ________________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_la.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_la.py:4: in
from shenfun.chebyshev.la import PDMA
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
_______________ ERROR collecting tests/test_lagrangian_particles.py _______________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_lagrangian_particles.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_lagrangian_particles.py:4: in
from shenfun import *
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
_____________________ ERROR collecting tests/test_matrices.py _____________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_matrices.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_matrices.py:9: in
import shenfun
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
________________ ERROR collecting tests/test_tensorproductspace.py ________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_tensorproductspace.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_tensorproductspace.py:8: in
from shenfun.fourier import bases as fbases
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
____________________ ERROR collecting tests/test_transforms.py ____________________
ImportError while importing test module '/home/fpoulin/software/shenfun/tests/test_transforms.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
tests/test_transforms.py:7: in
import shenfun
shenfun/init.py:33: in
from . import chebyshev
shenfun/chebyshev/init.py:2: in
from . import la
shenfun/chebyshev/la.py:6: in
from shenfun.optimization.cython import la
shenfun/optimization/cython/init.py:1: in
from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve,
E ModuleNotFoundError: No module named 'shenfun.optimization.cython.la'
!!!!!!!!!!!!!!!!!!!! Interrupted: 10 errors during collection !!!!!!!!!!!!!!!!!!!!!
=============================== 10 errors in 0.66s ================================`

parallel netcdf4

It'd be great to write netcdf output files.
This is the standard for many users in the geophysical community.
Parallel file writing is implemented in netCDF4.
I failed to install it via conda however.
Could it be possible to have it in the spectralDNS channel as it is the case for hdf5-parallel?

computational performance

I have not seen either in the shenfun paper or on spectralDNS/shenfun
guidelines regarding large scale applications, e.g: what operations are costly at the end?
what variables hold data at the end?
Note that I realize your paper may have what is required to answer these questions.
Yet after reading the paper, the answers are not straightforward for me.

It's maybe just an issue of form, here are some suggestions about form:

a cheat sheet with dos and don'ts related to code writing/architecture/variable handlings
code samples along with their profiling (can we use snakeviz in hpc configurations? that would be cool)
a feeback from your experience with your DNS setups

In my case for example, I am not familiar with spectral methods and their bottlenecks.
I am also always nervous about hidden python overheads or subtleties that may impact performance.
At the stage I'm at, i.e. planning the development of code based on shenfun that will ultimately run
at hpc scale, this seems important.

SymPy Deprecation

Hi,
very nice platform for spectral methods.
I have just installed following the instructions for conda and I found that Sympy gives a warning from

shenfun/lib/python3.9/site-packages/sympy/utilities/lambdify.py:807: SymPyDeprecationWarning:

The list of arguments is a set. This leads to unpredictable results
has been deprecated since SymPy 1.6.3. Use : Convert set into list or
tuple instead. See sympy/sympy#20013 for
more info.

for some demos (dirichlet_poisson2D, dirichlet_poisson3D, dirichlet_Helmholtz2D), but I did all the demos.

Best,
Alberto

import shenfun: Illegal instruction (core dumped) error

Dear mikaem,

I am trying to install shenfun package through anaconda 3.5.0 with python-3.6.0
but it is not working and giving core dumped on exporting it,

please help me to resolve the issue.

$ python
Python 3.6.0 | packaged by conda-forge | (default, Feb 9 2017, 14:36:55)
[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux
Type "help", "copyright", "credits" or "license" for more information.

import shenfun
Illegal instruction (core dumped)

shallow water 1 layer, issues with gradient

Hi,

I am trying to build upon KleinGordon in order to setup a 1 layer shallow water model:
https://github.com/apatlpo/shenfun_trial/blob/master/swater_1L/swater_1L.py

Unfortunately, the time rate of change of variables is very small at the moment and I suspect I am not normalizing correctly terms corresponding to g\partial_x h.
In other terms, I am missing the A term in the KleinGordon example.

If anybody has an idea, I would be most grateful.

Chebysehv and Fourier not working after installation

Hi,

I have started using shenfun. I installed the code using pip. I am trying to test some of the code given in the demo section. It seems to me that the Fourier and Chebyshev not working after the installation. Following things, I did on the Linux(Ubuntu) terminal. I am using python-3.6.

from shenfun import Basis
F0 = Basis(16, 'F')
Traceback (most recent call last):
File "", line 1, in
File "/home/aditya0006/.local/lib/python3.6/site-packages/shenfun/forms/arguments.py", line 84, in Basis
return B(N, **par)
File "/home/aditya0006/.local/lib/python3.6/site-packages/shenfun/fourier/bases.py", line 332, in init
self.plan((int(padding_factor*N),), (0,), np.float, {})
File "/home/aditya0006/.local/lib/python3.6/site-packages/shenfun/fourier/bases.py", line 277, in plan
xfftn_fwd = plan_fwd(U, s=s, axes=axis, threads=threads, flags=flags)
File "/home/aditya0006/.local/lib/python3.6/site-packages/mpi4py_fft/fftw/xfftn.py", line 240, in rfftn
flags, 1.0/M)
File "/home/aditya0006/.local/lib/python3.6/site-packages/mpi4py_fft/fftw/factory.py", line 102, in get_planned_FFT
assert dtype.upper() in fftlib
AssertionError

I tried a fresh install. The error is repeating.

Any help would be really great.

Thanks
Subhajit

issues with demo/KleinGordon.py out of the box

These are really minor issues but may be of some help.

I installed shenfun with:

conda create --name shenfun -c conda-forge -c spectralDNS python=3.6 shenfun

h5py was not installed and lead to a first crash while launching mpirun -np 4 python KleinGordon.py
after installing h5py (with conda), a second crash occured:

  File "KleinGordon.py", line 25, in <module>
    from spectralDNS.utilities import Timer
ModuleNotFoundError: No module named 'spectralDNS'

after commenting out the import statement and the call to Timer(), I got another crash:

Traceback (most recent call last):
  File "KleinGordon.py", line 59, in <module>
    file0 = HDF5Writer("KleinGordon{}.h5".format(N[0]), ['u', 'f'], TT)
  File "/Users/aponte/.miniconda3/envs/shenfun/lib/python3.6/site-packages/shenfun/utilities/h5py_writer.py", line 25, in __init__
    self.f = h5py.File(h5name, "w", driver="mpio", comm=T.comm)
  File "/Users/aponte/.miniconda3/envs/shenfun/lib/python3.6/site-packages/h5py/_hl/files.py", line 268, in __init__
    fapl = make_fapl(driver, libver, **kwds)
  File "/Users/aponte/.miniconda3/envs/shenfun/lib/python3.6/site-packages/h5py/_hl/files.py", line 71, in make_fapl
    kwds.setdefault('info', mpi4py.MPI.Info())
NameError: name 'mpi4py' is not defined

This one looks strange to me, mpi4run was installed with shenfun.
A quick search online seems to indicate h5py install with conda may be faulty.
Any suggestions?

Assembling matrix failed

Hi Mikael,

I am solving a problem of linear gradient elasticity in 2D. When I try to assemble one of the matrices I get a KeyError and I'm not sure what's supposed to be the problem. The code below should explain the problem:

from sympy import symbols, cos, sin, pi
from mpi4py import MPI
from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, Array,
        FunctionSpace, TensorProductSpace, VectorTensorProductSpace
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD
family = 'legendre'
x, y = symbols("x,y")

# computational domain
l = 140
d = l/2
domain_x = (0, l)
domain_y = (0, d)

# Size of discretization
n = 40
N = [n, n]

# spaces
BX_u_x = FunctionSpace(N[0], family=family, domain=domain_x, bc=(0, 1, 0, 0))
BY_u_x = FunctionSpace(N[1], family=family, domain=domain_y, bc=None)
BX_u_y = FunctionSpace(N[0], family=family, domain=domain_x, bc=(0, 0, 0, 0))
BY_u_y = FunctionSpace(N[1], family=family, domain=domain_y, bc=None)

subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
TX = TensorProductSpace(subcomms, (BX_u_x, BY_u_x), axes=(1, 0))
TY = TensorProductSpace(subcomms, (BX_u_y, BY_u_y), axes=(1, 0))
V = VectorTensorProductSpace([TX, TY])

u = TrialFunction(V)
v = TestFunction(V)

# assembling this matrix doesn't work
A = inner(div(grad(u)), div(grad(v)))

The error I get is the following one:

KeyError: -2

Best regards,

Thilo

An error related to the MPI communicators

Hello,

Thanks for this great library.

I have a question. I have created a solver with shenfun in a function which is used in another code in an iterative process. However, it only works with a limited number of iterations, and if the number of iterations goes up, the shenfun gives the following error.

mpi4py/MPI/Comm.pyx in mpi4py.MPI.Cartcomm.Sub()

Exception: Other MPI error, error stack:
PMPI_Cart_sub(213)..................: MPI_Cart_sub(comm=0x84000000, remain_dims=0x7f1aedab1320, comm_new=0x7f1aed8c7da0) failed
PMPI_Cart_sub(152)..................: 
MPIR_Comm_split_impl(253)...........: 
MPIR_Get_contextid_sparse_group(602): Too many communicators (0/2048 free on this process; ignore_id=0)

I think I should do something like unload the MPI or free the MPI communicators after each iteration (each time solver is used). Is there any way to solve this problem?

Sincerely

conda install mpi4py-fft

This may be a trivial question but: why is it that, when I install shenfun via conda, mpi4py-fft is not installed even though it is listed as a dependency in the install doc ?

Is it possible to have full periodic boundaris?

I'm trying to solve 2D Poisson equation with the periodic boundaries for both X and Y directions. I took the code from 3D example with two periodic dimensions and modified it accordingly, but it gives me an an error, which I can't understand.

The code is the following:

import sys, os
import importlib
from sympy import symbols, cos, sin, lambdify
import numpy as np
from shenfun import inner, div, grad, TestFunction, TrialFunction, Array, \
    Function, Basis, TensorProductSpace
import time
from mpi4py import MPI


comm = MPI.COMM_WORLD

base = importlib.import_module('.'.join(('shenfun', 'legendre')))
Solver = base.la.Helmholtz

data = np.loadtxt('q_distrib.dat')
fe = np.array(data)

print(np.sum(fe))
fe -= np.mean(fe)
print(np.sum(fe))

# Size of discretization
N = fe.shape
print('N=',N)

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

# Get f on quad points
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)

# Get left hand side of Poisson equation
matrices = inner(v, div(grad(u)))

# Create Helmholtz linear algebra solver
H = Solver(*matrices)

# Solve and transform to real space
u_hat = Function(T)           # Solution spectral space
u_hat = H(u_hat, f_hat)       # Solve
uq = u_hat.backward()
uh = uq.forward()

np.savetxt('solution_2d.dat',uq)

The following error is produced:

Traceback (most recent call last):
  File "shen2d.py", line 29, in <module>
    T = TensorProductSpace(comm, (K1, K2), axes=(0,1))
  File "/home/syesylev/.local/lib/python3.6/site-packages/shenfun/tensorproductspace.py", line 158, in __init__
    xfftn.plan(pencilB.subshape, axes, dtype, kw)
  File "/home/syesylev/.local/lib/python3.6/site-packages/shenfun/fourier/bases.py", line 277, in plan
    xfftn_fwd = plan_fwd(U, s=s, axes=axis, threads=threads, flags=flags)
  File "/home/syesylev/.local/lib/python3.6/site-packages/mpi4py_fft/fftw/xfftn.py", line 229, in rfftn
    assert input_array.dtype.char in 'fdg'
AssertionError

Apparently I'm missing something but I can't figure out what exactly.

Convergence troubles

Hi Mikael,

I have been contacting you about half a year ago regarding the use of shenfun for solving problems in linear elasticity and linear gradient elasticity. I am facing some troubles while performing a convergence test with a manufactured solution. To keep it simple, I am only describing the case of standard linear elasticity. I would greatly appreciate it if you could take a look at the code for the problem described below and give us your opinion.

I have derived two different variational formulations for the spectral Galerkin method. The first one is based on the principle of virtual work which reads

$\int_{\Omega} \boldsymbol{\sigma} : grad(\delta \boldsymbol{v}) d V = \int_{\Omega} \rho \boldsymbol{b} \cdot \delta \boldsymbol{v} d V + \int_{\partial \Omega} \boldsymbol{t} \cdot \delta \boldsymbol{v} d A \,.$
Inserting Hooke's Law

$\boldsymbol{\sigma} = \lambda \nabla \cdot \boldsymbol{u} \boldsymbol{1} + \mu (\nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^T)$
yields the variational form

$\int_\Omega ( \lambda (\nabla \cdot \boldsymbol{u}) ( \nabla \cdot \delta \boldsymbol{v}) + \frac{\mu}{2} (\nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^T) : \nabla \delta \boldsymbol{v}) d V$

$= \int \limits_{\Omega} \rho \boldsymbol{b} \cdot \delta \boldsymbol{v} d V + \int_{\partial \Omega} \boldsymbol{t} \cdot \delta \boldsymbol{v} d A$

My manufactured solution only contains Dirichlet boundary conditions. Therefore, the test functions vanish on entire boundary and so the boundary integral.

The second variational formulation is based on the Lamé-Navier equations which reads

$(\lambda + \mu) \nabla (\nabla \cdot \boldsymbol{u}) + \mu \Delta \boldsymbol{u} + \rho \boldsymbol{b} = \boldsymbol{0}$
Through the integration by part the weak form is obtained as

$\int_{\Omega} (\lambda + \mu) (\nabla \cdot \boldsymbol{u}) (\nabla \cdot \delta \boldsymbol{v}) + \mu (\nabla \boldsymbol{u} ) : (\nabla \delta \boldsymbol{v}) d V$

$= \int_{\Omega} \rho \boldsymbol{b} \cdot \delta \boldsymbol{v} d V + \int_{\partial \Omega} \delta \boldsymbol{v} \cdot ( \mu \nabla \boldsymbol{u} + (\lambda + \mu) (\nabla \cdot \boldsymbol{u}) \boldsymbol{1}) \cdot \boldsymbol{n} d A \,.$
For the manufactured solution (pure Dirichlet problem) the boundary terms again vanishes.

I think that both variational formulations are equivalent for a pure Dirichlet problem. In order to check that both formulations are working in shenfun, I performed a convergence test. I computed the integral error of the numerical w.r.t. the manufactured one. I found that the error goes down to 1e-15 for N=30 basis functions in case of the second version. For the principle of virtual power (first version), however, the error is still around 1e-3.
is working where as the first one does not. Maybe you could take a look at the attached code file and give me your opinion what goes wrong and how it could possibly be fixed.

The reason for considering the first version at all is that I want to study another problem which is not a pure Dirichlet problem. Instead it contains one homogeneous traction boundary condition. In this problem, the boundary integral above also vanishes but the second variational formulation is not correct since the boundary integral does not contain the (physically correct) traction.

The code below should explain the problem. Both variational forms are inplemented in the code.

from sympy import symbols, cos, sin, pi
from math import sqrt
from shenfun import comm, TensorProductSpace, VectorSpace, FunctionSpace, Function, inner, Array, TrialFunction, \
    TestFunction, grad, div, Dx, BlockMatrix, extract_bc_matrices, project

# computational domain
l = 1.
h = l/2
domain_x = (0, l)
domain_y = (0, h)
domain = (domain_x, domain_y)
dim = len(domain)
# size of discretization
N = 30

# elastic constants
E = 400. # Young's modulus
nu = 0.4 # Poisson's ratio
# Lame parameters
lambd = E * nu / ( (1.0 + nu) * (1.0 - 2.0*nu) )
mu = E / (2.0 * (1.0 + nu) )

# displacement value
u0 = 1.

# sympy variables
x, y = symbols("x,y")
coord = [x, y]

# manufactured solution (sympy expression)
ua_sp = (u0*((1+x/l)*(y/h)**2*(1-y/h)**2*sin(2*pi*x/l)*cos(3*pi*y/h) + x/l*4*y/h*(1-y/h)), u0*x/l*(1-x/l)*sin(2*pi*y/h))

# Section 1) Setup  boundary conditions, shenfun function spaces etc.
# 1.1) Specify boundary conditions for the displacement
bcs = (
       # first component
       (
        # (left, right)
        (0., u0*4.0 * y / h * (1.0 - y/h) ),
        # (bottom, top)
        (0., 0.)
       ),
       # second component
       (
        # (left, right)
        (0., 0.),
        # (bottom, top)
        (0., 0.)
       )
      )

# 1.2) Create a VectorSpace for the displacement
vector_tensor_spaces = []
for i in range(dim): # nb of displacement components
    tensor_spaces = []
    for j in range(dim): # nb of FunctionSpaces for each component
        basis = FunctionSpace(N, family='legendre', bc=bcs[i][j], domain=domain[j])
        tensor_spaces.append(basis)
    vector_tensor_spaces.append(TensorProductSpace(comm, tuple(tensor_spaces)))
V = VectorSpace(vector_tensor_spaces)

# 1.3) Test and trial functions
u = TrialFunction(V)
v = TestFunction(V)

# 1.4) Compute the body force
# create a VectorSpace without boundary conditions
V_none = V.get_orthogonal()

# compute analytical solution as shenfun Function
ua = Array(V_none, buffer=ua_sp).forward()

# compute symbolic body force (as sympy function)
sym_body_force = (-(lambd+mu)*grad(div(ua))-mu*div(grad(ua))).tosympy(basis=ua_sp, psi=(x, y))

# compute body forces at quadrature points
body_force = Array(V_none, buffer=tuple(sym_body_force))

# Section 2) Solve the problem

# 2.1) Define weak form and matrices
matrices = []

#------------------------------------------------------------------------------#
# 1st version: principle of virtual power
#------------------------------------------------------------------------------#
A = inner(mu*grad(u), grad(v))

# inner(grad(u)^T, grad(v))
B = []
for i in range(dim):
    for j in range(dim):
            temp = inner(mu*Dx(u[i], j), Dx(v[j], i))
            if isinstance(temp, list):
                B += temp
            else:
                B += [temp]

C = inner(lambd*div(u), div(v))

matrices = A + B + C

#------------------------------------------------------------------------------#
# 2nd version: Lame-Navier equations
#------------------------------------------------------------------------------#
#A = inner(mu*grad(u), grad(v))
#
#B = inner((lambd + mu)*div(u), div(v))
#
#matrices = A + B

# 2.2) Right-hand side from body forces
b = inner(v, body_force)

# solution
u_hat = Function(V)

# 2.3) Setup the linear system 
# get boundary matrices
bc_mats = extract_bc_matrices([matrices])

# BlockMatrix for homogeneous part
M = BlockMatrix(matrices)

# BlockMatrix for inhomogeneous part
BM = BlockMatrix(bc_mats)

# inhomogeneous part of solution
uh_hat = Function(V).set_boundary_dofs()        
# additional part to be passed to the right hand side
b_add = Function(V)
b_add = BM.matvec(-uh_hat, b_add) # negative because added to right hand side

# 2.4) Solve the linear system 
# solve homogeneous part of solution
u_hat = M.solve(b + b_add)

# solution
u_hat += uh_hat

# Section 3) Compute error

# compute integral error
error_array = Array(V, buffer=ua_sp)
error_array -= project(u_hat, V).backward()
error_ua = sqrt(inner((1, 1), error_array**2))

print('error (analytical solution):', error_ua)

Best regards,

Thilo

2D example crashes with error

When trying to run the dirichlet_poisson2D.py the following error appears:

$ python3 dirichlet_poisson2D.py 64 legendre
Traceback (most recent call last):
  File "dirichlet_poisson2D.py", line 49, in <module>
    K1 = Basis(N[1], family='F', dtype='d', domain=(-2*np.pi, 2*np.pi))
  File "/home/syesylev/.local/lib/python3.6/site-packages/shenfun/forms/arguments.py", line 84, in Basis
    return B(N, **par)
  File "/home/syesylev/.local/lib/python3.6/site-packages/shenfun/fourier/bases.py", line 332, in __init__
    self.plan((int(padding_factor*N),), (0,), np.float, {})
  File "/home/syesylev/.local/lib/python3.6/site-packages/shenfun/fourier/bases.py", line 277, in plan
    xfftn_fwd = plan_fwd(U, s=s, axes=axis, threads=threads, flags=flags)
  File "/home/syesylev/.local/lib/python3.6/site-packages/mpi4py_fft/fftw/xfftn.py", line 240, in rfftn
    flags, 1.0/M)
  File "/home/syesylev/.local/lib/python3.6/site-packages/mpi4py_fft/fftw/factory.py", line 102, in get_planned_FFT
    assert dtype.upper() in fftlib
AssertionError

An example file is not modified what is the problem?

warnings on test_io

This is not a bug but when I test test_io I get a warning, which I presume you would want to fix up sometime?

`(shenfun) fpoulin@amlt21:~$ python -m pytest software/shenfun/tests/test_io.py
=============================== test session starts ===============================
platform linux -- Python 3.8.1, pytest-5.3.5, py-1.8.1, pluggy-0.13.1
rootdir: /home/fpoulin/software/shenfun
collected 24 items

software/shenfun/tests/test_io.py ........................ [100%]

================================ warnings summary =================================
software/anaconda3/envs/shenfun/lib/python3.8/site-packages/shenfun/forms/arguments.py:700
/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/shenfun/forms/arguments.py:700: PytestCollectionWarning: cannot collect test class 'TestFunction' because it has a init constructor (from: tests/test_io.py)
class TestFunction(BasisFunction):

tests/test_io.py::test_regular_2D[hdf5-False]
tests/test_io.py::test_mixed_2D[True-hdf5-False]
tests/test_io.py::test_mixed_2D[False-hdf5-False]
tests/test_io.py::test_regular_3D[hdf5-False]
tests/test_io.py::test_mixed_3D[True-hdf5-False]
tests/test_io.py::test_mixed_3D[False-hdf5-False]
/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/mpi4py_fft/io/generate_xdmf.py:117: H5pyDeprecationWarning: The default file mode will change to 'r' (read-only) in h5py 3.0. To suppress this warning, pass the mode you need to h5py.File(), or set the global default h5.get_config().default_file_mode, or set the environment variable H5PY_DEFAULT_READONLY=1. Available modes are: 'r', 'r+', 'w', 'w-'/'x', 'a'. See the docs for details.
f = h5py.File(h5filename)

-- Docs: https://docs.pytest.org/en/latest/warnings.html
========================= 24 passed, 7 warnings in 1.35s ==========================
`

controlling numerical instabilities and masks

I have developed a 2D code to solve the vorticity equation using a pseudo-spectral approach. It works beautifully in the case when it is Fourier by Fourier, if i use an exponential filter to tame the numerical instability. I'm now working on a cheb by Fourier version, and realizing I have a lot to learn on how shenfun works and I am hoping to get some guidance.

My first thought was to essentially truncate the highest Fourier and cheb modes. I was hoping that mask might do this but I see that when I create the mask for a Function space, it gives me a 1D vector for the mask in the Fourier direction. This would only control part of the issue and I need to do something to control the small scale features in both directions. If I were to build a 2D mask function that was 1 for the first 2/3 of the domain in each direction, and multiply the transformed flux by this masked function, do you think that would do the trick?

I also looked into padding but that seems to be most useful when multiplying functions and not necessarily for filtering the resulting flux, as far as I can tell.

Any advice would be greatly apprecitaed.

Reference for IRK3 integrator

Hi,

I was just curious about what particular kind of RK scheme was being used in the IRK3 integrator, and if there was a particular reference from which you drew the coefficients?

Thanks.

test_io.py error; mpi version problem?

Hi,

When I run:
bash rundemos.sh

all the demos run successfully. However, when I run:
python -m pytest .

(from shenfun/tests/ folder)

I run into the MPI error (attached below).

I used:
conda create --name shenfun_prec -c conda-forge -c spectralDNS shenfun

to create this shenfun environment in conda.

I also tried building the environment using instructions here:
https://shenfun.readthedocs.io/en/latest/installation.html

and I run into the same exact error.

Moreoever, running (after installing spectralDNS):
mpirun -np 4 python TG.py NS

for spectralDNS ends up with the same error as well (spectralDNS library, not shenfun).

Is it a way to fix this?

Thanks,
Farrukh

============================= test session starts ==============================
platform linux -- Python 3.6.7, pytest-4.1.1, py-1.7.0, pluggy-0.8.1
rootdir: /home/fnauman/python3/shenfun, inifile:
collected 3539 items / 1 errors

==================================== ERRORS ====================================
______________________ ERROR collecting tests/test_io.py _______________________
test_io.py:13: in
import h5py
../../../anaconda3/envs/shenfun_prec/lib/python3.6/site-packages/h5py/init.py:24: in
from . import _errors
MPI.pxd:28: in init h5py._errors
???
E ValueError: mpi4py.MPI.Status has the wrong size, try recompiling. Expected 48, got 40
=============================== warnings summary ===============================
/home/fnauman/anaconda3/envs/shenfun_prec/lib/python3.6/site-packages/shenfun/forms/arguments.py:510
/home/fnauman/anaconda3/envs/shenfun_prec/lib/python3.6/site-packages/shenfun/forms/arguments.py:510: PytestWarning: cannot collect test class 'TestFunction' because it has a init constructor
class TestFunction(BasisFunction):

-- Docs: https://docs.pytest.org/en/latest/warnings.html
!!!!!!!!!!!!!!!!!!! Interrupted: 1 errors during collection !!!!!!!!!!!!!!!!!!!!
===================== 1 warnings, 1 error in 1.74 seconds ======================

Chebyshev Biharmonic BCs

Hi,

very much like this Project!

Here is a small bug I've encountered
Using biharmonic BCs in combination with Chebyshev bases (like in the Rayleigh Benard demo), shenfun raises a not implemented error:

from shenfun import *
VB = FunctionSpace(100, 'Chebyshev', bc="biharmonic")

I fixed it by adding string recognition to the chebyshev bases part in forms/arguments.py (line 234):

    if isinstance(bc, (tuple, dict)):
        key, par['bc'] = _process_bcs(bc, domain)

    elif bc is None:
        key = ''
    #------------------------------------------
    elif isinstance(bc, str):
        if bc.lower() == 'biharmonic':  
            key = "LDLNRDRN"
        else:
            raise NotImplementedError
    #------------------------------------------
    else:
        raise NotImplementedError

Best regards

ImportError: SolverGeneric1ND

Hello Mikael,

I hope you are well.

After not using shenfun for a while I picked it picked it up today and see you've done some really interesting things, like including some demos in polar coordinates. They look great!

Unfortunately, when I try the new code, for example unitdisk_helmholtz.py, I get a problem with importing Solvergeneric1ND.

I did a conda update and thought that would help but it didn't. Do you recommend that I do a fresh install and that should add it? Or is there an additional step that I need to do?

Thanks in advance,
Francis

2D Rayleigh-Benard in vorticity-streamfunction formulation

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 :

$\begin{aligned} \frac{\partial \omega}{\partial t} + J(\psi, \omega) & = -Ra Pr \frac{\partial \theta}{\partial x} + Pr \nabla^2 \omega \\ \frac{\partial \theta}{\partial t} + J(\psi, \theta) & = \frac{\partial \psi}{\partial x} + \nabla^2 \theta \\ -\nabla^2 \psi & = \omega \end{aligned}$

where $\omega$ is the out-of-plane vorticity, $\psi$ is the streamfunction and $\theta$ 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 $\theta$ on both walls
Bihamonic boundary conditions for $\psi$ 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 $Q = \begin{bmatrix} \omega & \theta \end{bmatrix}^T$ 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

$\mathcal{L} = \begin{bmatrix} Pr \nabla^2 & 0 \\ 0 & \nabla^2 \end{bmatrix}$

( where I assume that I treat $-Ra Pr \frac{\partial \theta}{\partial x}$ and $\frac{\partial \psi}{\partial x}$ 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.

TypeError when passing a lambdified function as a buffer

Hi Mikael,

I am reporting another issue which is related to the evaluation of a scalar function on a 2D grid at the quadrature points. I used the same code as in the 3D Poisson's equation demo. However, in the case that the function doesn't depend on the entire set of spatial variables, I get a TypeError if I am using the lambdified version. I found a way to circumvent the error by passing the symbolic expression. I hope that the code below shows what I mean. However, I want to let you know that the problem might cause some headaches for inexperienced users like me.

Best regards,
Thilo

from sympy import symbols, lambdify
from mpi4py import MPI
from shenfun import Array, Basis, TensorProductSpace
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD
family = 'legendre'

# Use sympy
x, y = symbols("x,y")
# symbolic expression
fe_working = x*(x+y)
fe_non_working = x
# lambdify
fl_working = lambdify((x, y), fe_working, 'numpy')
fl_non_working = lambdify((x, y), fe_non_working, 'numpy')

# Size of discretization
n = 25
N = [n, n+1]

BX = Basis(N[0], family=family, bc=(0, 0))
BY = Basis(N[1], family=family, bc=(0, 0))
subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T = TensorProductSpace(subcomms, (BX, BY), axes=(1, 0))
X = T.local_mesh()

# Get f on quad points

# working case with expression
fj = Array(T, buffer=fe_working)
fj = Array(T, buffer=fe_non_working)

# working case with lambdify
fj = Array(T, buffer=fl_working(*X))

# non working case with lambdify
fj = Array(T, buffer=fl_non_working(*X)) # TypeError: buffer is too small for requested array

testing parallelization

I have tried to test my code to see how efficient the parallelization is in going from 1,2,4,8 and 16 cores and I've found some odd behaviour. Even if I use mpiexec -np 1, top says that it can use 4 or 5 cores. In general this would be a bonus, but it means that I can't really test the scalability.

Why is it that my shenfun, or more generally, python script can use more than the number of cores I specify?

Is there a way to enforce this request?

Failed to install on Linux Mint

I used shenfun before with great success but I can't install it on my new machine running pretty standard Linux Mint distro.
Below is the error message.
What should I do?

 pip3 install shenfun
Collecting shenfun
  Using cached shenfun-3.1.2.tar.gz (174 kB)
Requirement already satisfied: cython in /usr/lib/python3/dist-packages (from shenfun) (0.29.14)
Requirement already satisfied: mpi4py in /usr/lib/python3/dist-packages (from shenfun) (3.0.3)
Requirement already satisfied: mpi4py-fft in /home/semen/.local/lib/python3.8/site-packages (from shenfun) (2.0.3)
Requirement already satisfied: numpy in /usr/lib/python3/dist-packages (from shenfun) (1.17.4)
Requirement already satisfied: scipy in /usr/lib/python3/dist-packages (from shenfun) (1.3.3)
Building wheels for collected packages: shenfun
  Building wheel for shenfun (setup.py) ... error
  ERROR: Command errored out with exit status 1:
   command: /usr/bin/python3 -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"'; __file__='"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' bdist_wheel -d /tmp/pip-wheel-cwo1m0fk
       cwd: /tmp/pip-install-zotpvdlq/shenfun/
  Complete output (175 lines):
  running bdist_wheel
  running build
  running build_py
  creating build
  creating build/lib.linux-x86_64-3.8
  creating build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/matrixbase.py -> build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/tensorproductspace.py -> build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/spectralbase.py -> build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/coordinates.py -> build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/la.py -> build/lib.linux-x86_64-3.8/shenfun
  copying shenfun/__init__.py -> build/lib.linux-x86_64-3.8/shenfun
  creating build/lib.linux-x86_64-3.8/shenfun/optimization
  copying shenfun/optimization/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization
  creating build/lib.linux-x86_64-3.8/shenfun/optimization/cython
  copying shenfun/optimization/cython/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/cython
  creating build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/biharmonic.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/helmholtz.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/tdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/pdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/chebyshev.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/fdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  copying shenfun/optimization/numba/legendre.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
  creating build/lib.linux-x86_64-3.8/shenfun/legendre
  copying shenfun/legendre/lobatto.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
  copying shenfun/legendre/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
  copying shenfun/legendre/la.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
  copying shenfun/legendre/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
  copying shenfun/legendre/bases.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
  package init file 'shenfun/legendre/fastgl/__init__.py' not found (or not a regular file)
  creating build/lib.linux-x86_64-3.8/shenfun/laguerre
  copying shenfun/laguerre/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
  copying shenfun/laguerre/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
  copying shenfun/laguerre/bases.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
  creating build/lib.linux-x86_64-3.8/shenfun/hermite
  copying shenfun/hermite/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
  copying shenfun/hermite/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
  copying shenfun/hermite/bases.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
  creating build/lib.linux-x86_64-3.8/shenfun/chebyshev
  copying shenfun/chebyshev/quasi.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
  copying shenfun/chebyshev/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
  copying shenfun/chebyshev/la.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
  copying shenfun/chebyshev/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
  copying shenfun/chebyshev/bases.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
  creating build/lib.linux-x86_64-3.8/shenfun/fourier
  copying shenfun/fourier/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
  copying shenfun/fourier/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
  copying shenfun/fourier/bases.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
  creating build/lib.linux-x86_64-3.8/shenfun/jacobi
  copying shenfun/jacobi/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
  copying shenfun/jacobi/la.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
  copying shenfun/jacobi/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
  copying shenfun/jacobi/bases.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
  creating build/lib.linux-x86_64-3.8/shenfun/forms
  copying shenfun/forms/arguments.py -> build/lib.linux-x86_64-3.8/shenfun/forms
  copying shenfun/forms/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/forms
  copying shenfun/forms/project.py -> build/lib.linux-x86_64-3.8/shenfun/forms
  copying shenfun/forms/operators.py -> build/lib.linux-x86_64-3.8/shenfun/forms
  copying shenfun/forms/inner.py -> build/lib.linux-x86_64-3.8/shenfun/forms
  creating build/lib.linux-x86_64-3.8/shenfun/utilities
  copying shenfun/utilities/lagrangian_particles.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
  copying shenfun/utilities/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
  copying shenfun/utilities/integrators.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
  creating build/lib.linux-x86_64-3.8/shenfun/io
  copying shenfun/io/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/io
  running build_ext
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.cpp
  warning: shenfun/optimization/cython/Matvec.pyx:1656:14: the result of using negative indices inside of code sections marked as 'wraparound=False' is undefined
  warning: shenfun/optimization/cython/Matvec.pyx:1656:22: the result of using negative indices inside of code sections marked as 'wraparound=False' is undefined
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.cpp
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.cpp
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp
  
  Error compiling Cython file:
  ------------------------------------------------------------
  ...
  # distutils: language = c++
  #cython: boundscheck=False
  #cython: wraparound=False
  #cython: language_level=3
  
  cimport fastgl_wrap
         ^
  ------------------------------------------------------------
  
  shenfun/legendre/fastgl/fastgl_wrap.pyx:6:8: 'fastgl_wrap.pxd' not found
  
  Error compiling Cython file:
  ------------------------------------------------------------
  ...
      N : int
          The total number of quadrature points
      k : int
          The k'th point in the N-point quadrature rule
      """
      f = fastgl_wrap.GLPair(N, k)
                    ^
  ------------------------------------------------------------
  
  shenfun/legendre/fastgl/fastgl_wrap.pyx:19:19: cimported module has no attribute 'GLPair'
  
  Error compiling Cython file:
  ------------------------------------------------------------
  ...
      N : int
          The total number of quadrature points
      k : int
          The k'th point in the N-point quadrature rule
      """
      f = fastgl_wrap.GLPair(N, k)
                    ^
  ------------------------------------------------------------
  
  shenfun/legendre/fastgl/fastgl_wrap.pyx:19:19: Compiler crash in AnalyseExpressionsTransform
  
  ModuleNode.body = StatListNode(fastgl_wrap.pyx:6:0)
  StatListNode.stats[2] = StatListNode(fastgl_wrap.pyx:9:0)
  StatListNode.stats[0] = DefNode(fastgl_wrap.pyx:9:0,
      doc = "Return point and weight k for N-point Gauss-Legendre rule\n\n    Parameters\n    ----------\n    N : int\n        The total number of quadrature points\n    k : int\n        The k'th point in the N-point quadrature rule\n    ",
      modifiers = [...]/0,
      name = 'getGLPair',
      np_args_idx = [...]/0,
      num_required_args = 2,
      outer_attrs = [...]/2,
      py_wrapper_required = True,
      reqd_kw_flags_cname = '0',
      used = True)
  File 'ExprNodes.py', line 5344, in infer_type: SimpleCallNode(fastgl_wrap.pyx:19:26,
      result_is_used = True,
      use_managed_ref = True)
  File 'ExprNodes.py', line 6832, in infer_type: AttributeNode(fastgl_wrap.pyx:19:19,
      attribute = 'GLPair',
      is_attribute = 1,
      needs_none_check = True,
      result_is_used = True,
      use_managed_ref = True)
  
  Compiler crash traceback from this point on:
    File "/usr/lib/python3/dist-packages/Cython/Compiler/ExprNodes.py", line 6832, in infer_type
      return node.entry.type
  AttributeError: 'NoneType' object has no attribute 'type'
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Cheb.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Cheb.c
  /usr/lib/python3/dist-packages/Cython/Compiler/Main.py:369: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Cheb.pyx
    tree = Parsing.p_module(s, pxd, full_module_name)
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/convolve.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/convolve.c
  /usr/lib/python3/dist-packages/Cython/Compiler/Main.py:369: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/convolve.pyx
    tree = Parsing.p_module(s, pxd, full_module_name)
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/outer.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/outer.c
  cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/applymask.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/applymask.c
  building 'shenfun.optimization.cython.Matvec' extension
  creating build/temp.linux-x86_64-3.8
  creating build/temp.linux-x86_64-3.8/tmp
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython
  x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.o -g0 -O3 -w -Ofast -ffast-math -march=native
  x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/Matvec.cpython-38-x86_64-linux-gnu.so -std=c++11
  building 'shenfun.optimization.cython.la' extension
  x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.o -g0 -O3 -w -Ofast -ffast-math -march=native
  x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/la.cpython-38-x86_64-linux-gnu.so -std=c++11
  building 'shenfun.optimization.cython.evaluate' extension
  x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.o -g0 -O3 -w -Ofast -ffast-math -march=native
  x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/evaluate.cpython-38-x86_64-linux-gnu.so -std=c++11
  building 'shenfun.legendre.fastgl.fastgl_wrap' extension
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre
  creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl
  x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.o -g0 -O3 -w -Ofast -ffast-math -march=native
  /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp:1:2: error: #error Do not use this file, it is the result of a failed Cython compilation.
      1 | #error Do not use this file, it is the result of a failed Cython compilation.
        |  ^~~~~
  error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
  ----------------------------------------
  ERROR: Failed building wheel for shenfun
  Running setup.py clean for shenfun
Failed to build shenfun
Installing collected packages: shenfun
    Running setup.py install for shenfun ... error
    ERROR: Command errored out with exit status 1:
     command: /usr/bin/python3 -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"'; __file__='"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-ptiffp84/install-record.txt --single-version-externally-managed --user --prefix= --compile --install-headers /home/semen/.local/include/python3.8/shenfun
         cwd: /tmp/pip-install-zotpvdlq/shenfun/
    Complete output (169 lines):
    running install
    running build
    running build_py
    creating build
    creating build/lib.linux-x86_64-3.8
    creating build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/matrixbase.py -> build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/tensorproductspace.py -> build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/spectralbase.py -> build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/coordinates.py -> build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/la.py -> build/lib.linux-x86_64-3.8/shenfun
    copying shenfun/__init__.py -> build/lib.linux-x86_64-3.8/shenfun
    creating build/lib.linux-x86_64-3.8/shenfun/optimization
    copying shenfun/optimization/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization
    creating build/lib.linux-x86_64-3.8/shenfun/optimization/cython
    copying shenfun/optimization/cython/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/cython
    creating build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/biharmonic.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/helmholtz.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/tdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/pdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/chebyshev.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/fdma.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    copying shenfun/optimization/numba/legendre.py -> build/lib.linux-x86_64-3.8/shenfun/optimization/numba
    creating build/lib.linux-x86_64-3.8/shenfun/legendre
    copying shenfun/legendre/lobatto.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
    copying shenfun/legendre/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
    copying shenfun/legendre/la.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
    copying shenfun/legendre/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
    copying shenfun/legendre/bases.py -> build/lib.linux-x86_64-3.8/shenfun/legendre
    package init file 'shenfun/legendre/fastgl/__init__.py' not found (or not a regular file)
    creating build/lib.linux-x86_64-3.8/shenfun/laguerre
    copying shenfun/laguerre/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
    copying shenfun/laguerre/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
    copying shenfun/laguerre/bases.py -> build/lib.linux-x86_64-3.8/shenfun/laguerre
    creating build/lib.linux-x86_64-3.8/shenfun/hermite
    copying shenfun/hermite/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
    copying shenfun/hermite/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
    copying shenfun/hermite/bases.py -> build/lib.linux-x86_64-3.8/shenfun/hermite
    creating build/lib.linux-x86_64-3.8/shenfun/chebyshev
    copying shenfun/chebyshev/quasi.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
    copying shenfun/chebyshev/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
    copying shenfun/chebyshev/la.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
    copying shenfun/chebyshev/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
    copying shenfun/chebyshev/bases.py -> build/lib.linux-x86_64-3.8/shenfun/chebyshev
    creating build/lib.linux-x86_64-3.8/shenfun/fourier
    copying shenfun/fourier/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
    copying shenfun/fourier/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
    copying shenfun/fourier/bases.py -> build/lib.linux-x86_64-3.8/shenfun/fourier
    creating build/lib.linux-x86_64-3.8/shenfun/jacobi
    copying shenfun/jacobi/matrices.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
    copying shenfun/jacobi/la.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
    copying shenfun/jacobi/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
    copying shenfun/jacobi/bases.py -> build/lib.linux-x86_64-3.8/shenfun/jacobi
    creating build/lib.linux-x86_64-3.8/shenfun/forms
    copying shenfun/forms/arguments.py -> build/lib.linux-x86_64-3.8/shenfun/forms
    copying shenfun/forms/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/forms
    copying shenfun/forms/project.py -> build/lib.linux-x86_64-3.8/shenfun/forms
    copying shenfun/forms/operators.py -> build/lib.linux-x86_64-3.8/shenfun/forms
    copying shenfun/forms/inner.py -> build/lib.linux-x86_64-3.8/shenfun/forms
    creating build/lib.linux-x86_64-3.8/shenfun/utilities
    copying shenfun/utilities/lagrangian_particles.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
    copying shenfun/utilities/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
    copying shenfun/utilities/integrators.py -> build/lib.linux-x86_64-3.8/shenfun/utilities
    creating build/lib.linux-x86_64-3.8/shenfun/io
    copying shenfun/io/__init__.py -> build/lib.linux-x86_64-3.8/shenfun/io
    running build_ext
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.cpp' Cython extension (up-to-date)
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.cpp' Cython extension (up-to-date)
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.cpp' Cython extension (up-to-date)
    cythoning /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.pyx to /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
    # distutils: language = c++
    #cython: boundscheck=False
    #cython: wraparound=False
    #cython: language_level=3
    
    cimport fastgl_wrap
           ^
    ------------------------------------------------------------
    
    shenfun/legendre/fastgl/fastgl_wrap.pyx:6:8: 'fastgl_wrap.pxd' not found
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
        N : int
            The total number of quadrature points
        k : int
            The k'th point in the N-point quadrature rule
        """
        f = fastgl_wrap.GLPair(N, k)
                      ^
    ------------------------------------------------------------
    
    shenfun/legendre/fastgl/fastgl_wrap.pyx:19:19: cimported module has no attribute 'GLPair'
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
        N : int
            The total number of quadrature points
        k : int
            The k'th point in the N-point quadrature rule
        """
        f = fastgl_wrap.GLPair(N, k)
                      ^
    ------------------------------------------------------------
    
    shenfun/legendre/fastgl/fastgl_wrap.pyx:19:19: Compiler crash in AnalyseExpressionsTransform
    
    ModuleNode.body = StatListNode(fastgl_wrap.pyx:6:0)
    StatListNode.stats[2] = StatListNode(fastgl_wrap.pyx:9:0)
    StatListNode.stats[0] = DefNode(fastgl_wrap.pyx:9:0,
        doc = "Return point and weight k for N-point Gauss-Legendre rule\n\n    Parameters\n    ----------\n    N : int\n        The total number of quadrature points\n    k : int\n        The k'th point in the N-point quadrature rule\n    ",
        modifiers = [...]/0,
        name = 'getGLPair',
        np_args_idx = [...]/0,
        num_required_args = 2,
        outer_attrs = [...]/2,
        py_wrapper_required = True,
        reqd_kw_flags_cname = '0',
        used = True)
    File 'ExprNodes.py', line 5344, in infer_type: SimpleCallNode(fastgl_wrap.pyx:19:26,
        result_is_used = True,
        use_managed_ref = True)
    File 'ExprNodes.py', line 6832, in infer_type: AttributeNode(fastgl_wrap.pyx:19:19,
        attribute = 'GLPair',
        is_attribute = 1,
        needs_none_check = True,
        result_is_used = True,
        use_managed_ref = True)
    
    Compiler crash traceback from this point on:
      File "/usr/lib/python3/dist-packages/Cython/Compiler/ExprNodes.py", line 6832, in infer_type
        return node.entry.type
    AttributeError: 'NoneType' object has no attribute 'type'
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Cheb.c' Cython extension (up-to-date)
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/convolve.c' Cython extension (up-to-date)
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/outer.c' Cython extension (up-to-date)
    skipping '/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/applymask.c' Cython extension (up-to-date)
    building 'shenfun.optimization.cython.Matvec' extension
    creating build/temp.linux-x86_64-3.8
    creating build/temp.linux-x86_64-3.8/tmp
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython
    x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.o -g0 -O3 -w -Ofast -ffast-math -march=native
    x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/Matvec.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/Matvec.cpython-38-x86_64-linux-gnu.so -std=c++11
    building 'shenfun.optimization.cython.la' extension
    x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.o -g0 -O3 -w -Ofast -ffast-math -march=native
    x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/la.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/la.cpython-38-x86_64-linux-gnu.so -std=c++11
    building 'shenfun.optimization.cython.evaluate' extension
    x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.o -g0 -O3 -w -Ofast -ffast-math -march=native
    x86_64-linux-gnu-g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/optimization/cython/evaluate.o -lm -o build/lib.linux-x86_64-3.8/shenfun/optimization/cython/evaluate.cpython-38-x86_64-linux-gnu.so -std=c++11
    building 'shenfun.legendre.fastgl.fastgl_wrap' extension
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre
    creating build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl
    x86_64-linux-gnu-gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/lib/python3/dist-packages/numpy/core/include -I/usr/include/python3.8 -c /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp -o build/temp.linux-x86_64-3.8/tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.o -g0 -O3 -w -Ofast -ffast-math -march=native
    /tmp/pip-install-zotpvdlq/shenfun/shenfun/legendre/fastgl/fastgl_wrap.cpp:1:2: error: #error Do not use this file, it is the result of a failed Cython compilation.
        1 | #error Do not use this file, it is the result of a failed Cython compilation.
          |  ^~~~~
    error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
    ----------------------------------------
ERROR: Command errored out with exit status 1: /usr/bin/python3 -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"'; __file__='"'"'/tmp/pip-install-zotpvdlq/shenfun/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-ptiffp84/install-record.txt --single-version-externally-managed --user --prefix= --compile --install-headers /home/semen/.local/include/python3.8/shenfun Check the logs for full command output.

Generalized Jacobi polynomials

Implement generalized Jacobi polynomials, specifically for applications to fractional PDEs.

mhd problem: floating precision issue?

I was trying to use MHD.py in (spectralDNS/solvers) with shenfun since that example uses shenfun to setup the problem anyway. I ran the Taylor Green problem and at the end there is an assertion error. The answer is turning out to be 0.125 while it should be:

        assert np.round(k - 0.124565408177, 7) == 0
        assert np.round(b - 0.124637762143, 7) == 0

The full file can be seen here:
https://gist.github.com/fnauman/e105bdf901616345b70145c6ded7c963#file-tgmhd_shenfun-py

On a related note: At the beginning of the NavierStokes.py demo, you mention

Not implemented for efficiency

How slow is the shenfun implementation compared to spectralDNS? (and why, if I may ask?)

Assembling matrix via inner(div(u), div(v)) failed

Hi Mikael,

I am trying to solve a linear elasticity problem on a 2D domain. I'm using Legendre polynomials and, thus, I have assembled the weak form using integration by parts, which contains to terms:

The inner product of the gradients of the test and trial functions inner(grad(u), grad(v))
The product of the divergences of the test and trial functions which I tried to implement using inner(div(u), div(v)).

The first term works fine. The second term inner(div(u), div(v)) raises an error. However, I found a workaround, see my code below. I think it would be convenient to be able to directly use inner(div(u), div(v)).

Best regards,
Thilo

from sympy import symbols, cos, sin
from mpi4py import MPI
from shenfun import inner, grad, div, Dx, TestFunction, TrialFunction, Array,\
    Basis, TensorProductSpace, VectorTensorProductSpace, BlockMatrix
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD

family = 'legendre'

# Size of discretization
n = 40
N = [n, n+1]

# Bases and spaces
BX_hom = Basis(N[0], family=family, bc=(0, 0))
BY_hom = Basis(N[1], family=family, bc=(0, 0))

subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T_hom = TensorProductSpace(subcomms, (BX_hom, BY_hom))

V = VectorTensorProductSpace([T_hom, T_hom])

# Function spaces
u = TrialFunction(V)
v = TestFunction(V)

# Compute matrices
# Assemble standard term
A = inner(grad(u), grad(v))

# Assemble (div(u), div(v))-term

# non-working
B = inner(div(u), div(v)) # raises an error

# workaround
C00 = inner(Dx(u[0], 0), Dx(v[0], 0))
C01 = inner(Dx(u[0], 0), Dx(v[1], 1))
C02 = inner(Dx(u[1], 1), Dx(v[0], 0))
C03 = inner(Dx(u[1], 1), Dx(v[1], 1))
C = [C00, C01, C02, C03]

# Compute BlockMatrix
M = BlockMatrix(A + C)

# Compute right hand side
# Use sympy to compute a rhs
x, y = symbols("x,y")
fe = (sin(x**2)*(1-y)*(1+y), y**2*(cos(x-y**2)))
fj = Array(V, buffer=fe)

b = inner(v, fj)

# Compute solution
u_hat = M.solve(b)
u_sol = u_hat.backward()

tests fail with cython optimization

I can happily say that all the tests pass without any optimization. Today I tried with cython and numba.

The good news is that with numba they all pass.

The bad news is that with cython I have 202 fail. This is both on my laptop and a computer cluster, but both are installed using conda so it's really the same system.

I will continue to use numba but wanted to point this out, in case you weren't aware of it already.

Incorrect solution of 1D Poisson equation

I'm struggling to understand why shenfun produces wrong magnitude of solution for simplest 1D Poisson equation. I can't figure out what I'm doing wrong, so any help is deeply appreciated.

Here is the code, which solved 1D Poissin equation first using shenfun and then by naive double integration:

import numpy as np
from shenfun import inner, div, grad, TestFunction, TrialFunction, \
    Array, Function, Basis


# Naive integration
def integrate(vec, dx):
    f = np.zeros(vec.shape)
    for i in range(1,len(vec)):
        f[i] = f[i-1] + (vec[i]+vec[i-1])*0.5*dx
    return f

#--------------------------------------
# Integrate with shenfun
#--------------------------------------

data = np.loadtxt('charge.xvg',comments=('#','@'))
fe = np.array(data[:,1])

# Size of discretization
N = fe.shape[0]

# Domain is set to what we have on X in data file
SD = Basis(N, family='legendre', bc=(0,0), domain=(data[0,0],data[-1,0]))
X = SD.mesh()

u = TrialFunction(SD)
v = TestFunction(SD)

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

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

# Get left hand side of Poisson equation
A = inner(v, div(grad(u)))

u_hat = Function(SD)
u_hat = A.solve(f_hat, u_hat)
uj = u_hat.backward()
uh = uj.forward()

with open('solution_shen.dat','w') as f:
    for x,v in zip(X,uj):
        f.write(f'{x} {v}\n')
    
#--------------------------------------
# Naive integration
#--------------------------------------
dx = data[1,0]-data[0,0] # step on X
f = integrate(data[:,1],dx) # Integrate once
u = integrate(f,dx)         # Integrade second time
with open('solution_naive.dat','w') as out:
    for i in range(len(u)):
        out.write(f'{data[i,0]} {u[i]}\n')

And here is the result:

I can't understand two things here:

Why magnitude of shenfun solution is obviously wrong. I'm using exactly the same data for both methods and the results should be identical.
Why the SD.mesh() returns different X points (see systematic shift of the curves on X) from what we have in data file despite correct number of points and domain?

Could somebody explain me what I'm missing and what should I modify to get the same result as naive integration?

defining basis functions + warning

I am learning shenfun and really impressed by it. However, I do have a couple of questions.

I first wrote a Fourier based solver for advection-diffusion and works very nicely. I've since tried to play with a Chebyshev version, since it allows for Dirichlet and Neumann boundary conditions. I have two questions.

First, is this how you suggest we define these functions in the latest version? I got this from old code and suspect I should update it but not sure if that is true.

# Defining basis functions F0 = ShenDirichletBasis (N[0], quad = 'GC', bc=(0.,0.), domain=(0.,L[0]), padding_factor=1, dealias_direct=True) F1 = ShenDirichletBasis (N[1], quad = 'GC', bc=(0.,0.), domain=(0.,L[1]), padding_factor=1, dealias_direct=True)

Second, I get the following warning.

`#error

 python Advection_shenfun.py 
/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site- 
packages/scipy/sparse/linalg/dsolve/linsolve.py:137: SparseEfficiencyWarning: spsolve requires A 
be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',

data layout

Hi,

I'm trying to understand the layout of data among MPI procs.
Could anybody please take a look at the following pics in order to see if I got it right?

Role of BCBasis and BoundaryValues

I would like to implement inhomogeneous mixed Dirichlet-Neumann BCs. Homogeneous mixed BCs are already available for example by using the class DirichletNeumannBasis. There are two types of classes which seem to handle the inhomogeneous Dirichlet BCs in shenfun: BCBasis and BoundaryValues. My question is whether both classes require a modification or if it is sufficient to modify BCBasis only.

Furthermore, it is not clear if both classes only work for the Dirichlet BCs or also for the mixed cases. All of the ansatz spaces are constructed such that the functions \phi_{N-2} and \phi_{N-1} fulfill the inhomogeneous BCs. All other functions fulfill homogeneous BCs. Thus, for u(x=-1)=a and u(x=1)=b, the coefficients are given by c_{N-2}=a and c_{N-1}=b. This relation regarding the coefficients also applies for the mixed BCs. Therefore, it seems to me that it is sufficient to modify BoundaryValues only.

I would be great if someone could clarify what is the purpose of the different classes such that I can try extend them for the mixed cases.

error with some demos ( + a question on diagnostics)

Hello Mikael,

I am looking at a couple of examples to figure out how to do diagnostics better than what I'm doing and I noticed that the two Ginzburg_Laudau examples don't work on my latest install. Also, when I tried biharmonic1D.py.

Below you will see that I did update shenfun and then you will see the three sets of errors that I have.

I will point out tha tthe tests all pass so I thought my system was good to go but given these errors maybe there is a problem.

After I can get these working I am curious to ask your advice on how best to have diagnostics in shenfun.

Cheers!

P.S. I see that I didn't have mayavi on my system so I tried to install it using pip3 but it failed to build. The essence seems to be in this statement here. Sadly, I still get an error tha tmayavi_show is not defined.

Failed to build mayavi
Installing collected packages: configobj, traits, pyface, traitsui, apptools, envisage, pygments, vtk, mayavi
    Running setup.py install for mayavi ... done
  DEPRECATION: mayavi was installed using the legacy 'setup.py install' method, because a wheel could not be built for it. pip 21.0 will remove support for this functionality. A possible replacement is to fix the wheel build issue reported above. You can find discussion regarding this at https://github.com/pypa/pip/issues/8368.

Downloading and Extracting Packages
pip-20.2.3           | 1.7 MB    | ############################################################################### | 100% 
setuptools-49.6.0    | 742 KB    | ############################################################################### | 100% 
attrs-20.2.0         | 42 KB     | ############################################################################### | 100% 
openssl-1.1.1h       | 2.5 MB    | ############################################################################### | 100% 
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

(shenfun) fpoulin@vortex2:~/software/shenfun/demo$ python Ginzburg_Landau.py 
No wisdom imported
Traceback (most recent call last):
  File "Ginzburg_Landau.py", line 91, in <module>
    U_hat = integrator.solve(U, U_hat, dt, (0, end_time))
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site-packages/shenfun/utilities/integrators.py", line 241, in solve
    self.U_hat0[:] = u_hat*self.ehL_h
ValueError: operands could not be broadcast together with shapes (129,129) (2,) 

(shenfun) fpoulin@vortex2:~/software/shenfun/demo$ python Ginzburg_Landau_sphere_IRK3.py 
Traceback (most recent call last):
  File "Ginzburg_Landau_sphere_IRK3.py", line 10, in <module>
    from mayavi import mlab
ModuleNotFoundError: No module named 'mayavi'
(

(shenfun) fpoulin@vortex2:~/software/shenfun/demo$ python biharmonic1D.py 24 chebyshev
Error=6.1129316770884646e-05
Traceback (most recent call last):
  File "biharmonic1D.py", line 88, in <module>
    assert np.linalg.norm(uj-uq) < 1e-8
AssertionError

matplotlib is missing in the conda package

Hi,
it occurred to me that matplotlib is missing in the conda package. It can of course be easily installed later. However it would be convenient if it is added to the package.
Best wishes,
Sebastian

software architecture around shenfun?

I will be in the next incoming months start writing a code library solving fluid flows with different sets of physics, in different types of domains (periodic vs closed).
The library will rely on shenfun and I was wondering whether anybody had advice regarding software architecture in order to do that.
One item I am worrying about for example is how to "gracefully" allow for switches between different spectral bases.

My strategy was to inspire myself from spectralDNS.
A minimal code library leveraging shenfun and providing an adequate architecture would actually be something useful to have around in order for people to start coding their solvers.

Does anybody has an opinion about this?

Stokes3D.py Incurs High Error after Changing the Coefficients of the Periodic Flow on the X-Y Plane

When one changes uex = sin(2*y)*(1-z**2) to uex = sin(2.05*y)*(1-z**2) in the manufactured solution, the error for the x-component of the flow velocity increases from 3.9229e-14 to 2.672e0.

Would I need to change some of the parameters when creating the Fourier basis for the periodic flow on the x-axis?

Furthermore, would it be possible to perform a singular value decomposition on the non-exact solution produced by Shenfun to determine which basis functions preserve the greatest variance of the fluid dynamics. Would the following code suffice?

r = 5
for dimension in range(0, 3):
            u_decomp, s_decomp, vh_decomp = np.linalg.svd(u_[dimension, :, :, :])
            print(u_decomp.shape, s_decomp.shape, vh_decomp.shape)
            for i in range(n):
                plt.figure()
                plt.yscale('log')
                plt.scatter(np.arange(n), s_decomp[i])
                plt.savefig("./sigma/" + str(dimension) + "_" + str(i) + "_variance_decomp.png")
            plt.figure()
            plt.yscale('linear')
            for i in range(0, r):
                plt.plot(np.arange(n), u_decomp[i])
            plt.savefig(str(dimension) + "_variance_modes.png")

variance_decomposition.png

variance_modes.png
The resulting modes are quite jumbled as I am having trouble separating the various axes of the block matrix M.

meshgrid.png
Here is my pictorial understanding of the local mesh X. However, I don't have an intuitive grasp on block matrix M.

Cython fails if shenfun is imported from conda

Hi,
I am trying to get the current github version of shenfun running in a conda-environment. I installed shenfun using conda as described in the documentation and cloned the github repository.

I could fix an error which was related to missing include files for numpy. But the include files for ios seem to be missing too, see below. I am somehow struggling to configure Cython such that the import shenfun works fine.

Best wishes,
Sebastian

(shenfun) sg@Piola:~/Documents/Projects/shenfun$ python
Python 3.8.2 | packaged by conda-forge | (default, Apr 24 2020, 08:20:52) 
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
>>> import pyximport
>>> pyximport.install(setup_args={"include_dirs":numpy.get_include()})
(None, <pyximport.pyximport.PyxImporter object at 0x7f854f0f95e0>)
>>> import shenfun
In file included from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/ndarraytypes.h:1832:0,
                 from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
                 from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /home/sg/.pyxbld/temp.linux-x86_64-3.8/pyrex/shenfun/optimization/cython/la.c:600:
/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it with " \
  ^~~~~~~
/home/sg/.pyxbld/temp.linux-x86_64-3.8/pyrex/shenfun/optimization/cython/la.c:602:10: fatal error: ios: No such file or directory
 #include "ios"
          ^~~~~
compilation terminated.
Traceback (most recent call last):
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 117, in _compile
    self.spawn(compiler_so + cc_args + [src, '-o', obj] +
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 910, in spawn
    spawn(cmd, dry_run=self.dry_run)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/spawn.py", line 36, in spawn
    _spawn_posix(cmd, search_path, dry_run=dry_run)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/spawn.py", line 157, in _spawn_posix
    raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 214, in load_module
    so_path = build_module(module_name, pyxfilename, pyxbuild_dir,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 186, in build_module
    so_path = pyxbuild.pyx_to_dll(pyxfilename, extension_mod,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyxbuild.py", line 102, in pyx_to_dll
    dist.run_commands()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 966, in run_commands
    self.run_command(cmd)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 985, in run_command
    cmd_obj.run()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 186, in run
    _build_ext.build_ext.run(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 340, in run
    self.build_extensions()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 195, in build_extensions
    _build_ext.build_ext.build_extensions(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 449, in build_extensions
    self._build_extensions_serial()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 474, in _build_extensions_serial
    self.build_extension(ext)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 528, in build_extension
    objects = self.compiler.compile(sources,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 574, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 120, in _compile
    raise CompileError(msg)
distutils.errors.CompileError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sg/Documents/Projects/shenfun/shenfun/__init__.py", line 33, in <module>
    from . import chebyshev
  File "/home/sg/Documents/Projects/shenfun/shenfun/chebyshev/__init__.py", line 2, in <module>
    from . import la
  File "/home/sg/Documents/Projects/shenfun/shenfun/chebyshev/la.py", line 6, in <module>
    from shenfun.optimization.cython import la
  File "/home/sg/Documents/Projects/shenfun/shenfun/optimization/cython/__init__.py", line 1, in <module>
    from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve, \
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 459, in load_module
    module = load_module(fullname, self.path,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 231, in load_module
    raise exc.with_traceback(tb)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 214, in load_module
    so_path = build_module(module_name, pyxfilename, pyxbuild_dir,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 186, in build_module
    so_path = pyxbuild.pyx_to_dll(pyxfilename, extension_mod,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyxbuild.py", line 102, in pyx_to_dll
    dist.run_commands()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 966, in run_commands
    self.run_command(cmd)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 985, in run_command
    cmd_obj.run()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 186, in run
    _build_ext.build_ext.run(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 340, in run
    self.build_extensions()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 195, in build_extensions
    _build_ext.build_ext.build_extensions(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 449, in build_extensions
    self._build_extensions_serial()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 474, in _build_extensions_serial
    self.build_extension(ext)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 528, in build_extension
    objects = self.compiler.compile(sources,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 574, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 120, in _compile
    raise CompileError(msg)
ImportError: Building module shenfun.optimization.cython.la failed: ["distutils.errors.CompileError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1\n"]
>>> import numpy
>>> import shenfun
In file included from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/ndarraytypes.h:1832:0,
                 from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
                 from /home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /home/sg/.pyxbld/temp.linux-x86_64-3.8/pyrex/shenfun/optimization/cython/la.c:600:
/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it with " \
  ^~~~~~~
/home/sg/.pyxbld/temp.linux-x86_64-3.8/pyrex/shenfun/optimization/cython/la.c:602:10: fatal error: ios: No such file or directory
 #include "ios"
          ^~~~~
compilation terminated.
Traceback (most recent call last):
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 117, in _compile
    self.spawn(compiler_so + cc_args + [src, '-o', obj] +
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 910, in spawn
    spawn(cmd, dry_run=self.dry_run)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/spawn.py", line 36, in spawn
    _spawn_posix(cmd, search_path, dry_run=dry_run)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/spawn.py", line 157, in _spawn_posix
    raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 214, in load_module
    so_path = build_module(module_name, pyxfilename, pyxbuild_dir,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 186, in build_module
    so_path = pyxbuild.pyx_to_dll(pyxfilename, extension_mod,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyxbuild.py", line 102, in pyx_to_dll
    dist.run_commands()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 966, in run_commands
    self.run_command(cmd)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 985, in run_command
    cmd_obj.run()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 186, in run
    _build_ext.build_ext.run(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 340, in run
    self.build_extensions()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 195, in build_extensions
    _build_ext.build_ext.build_extensions(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 449, in build_extensions
    self._build_extensions_serial()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 474, in _build_extensions_serial
    self.build_extension(ext)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 528, in build_extension
    objects = self.compiler.compile(sources,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 574, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 120, in _compile
    raise CompileError(msg)
distutils.errors.CompileError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sg/Documents/Projects/shenfun/shenfun/__init__.py", line 33, in <module>
    from . import chebyshev
  File "/home/sg/Documents/Projects/shenfun/shenfun/chebyshev/__init__.py", line 2, in <module>
    from . import la
  File "/home/sg/Documents/Projects/shenfun/shenfun/chebyshev/la.py", line 6, in <module>
    from shenfun.optimization.cython import la
  File "/home/sg/Documents/Projects/shenfun/shenfun/optimization/cython/__init__.py", line 1, in <module>
    from .la import TDMA_SymLU, TDMA_SymSolve, PDMA_SymLU, PDMA_SymSolve, \
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 459, in load_module
    module = load_module(fullname, self.path,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 231, in load_module
    raise exc.with_traceback(tb)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 214, in load_module
    so_path = build_module(module_name, pyxfilename, pyxbuild_dir,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyximport.py", line 186, in build_module
    so_path = pyxbuild.pyx_to_dll(pyxfilename, extension_mod,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/pyximport/pyxbuild.py", line 102, in pyx_to_dll
    dist.run_commands()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 966, in run_commands
    self.run_command(cmd)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/dist.py", line 985, in run_command
    cmd_obj.run()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 186, in run
    _build_ext.build_ext.run(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 340, in run
    self.build_extensions()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/site-packages/Cython/Distutils/old_build_ext.py", line 195, in build_extensions
    _build_ext.build_ext.build_extensions(self)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 449, in build_extensions
    self._build_extensions_serial()
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 474, in _build_extensions_serial
    self.build_extension(ext)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/command/build_ext.py", line 528, in build_extension
    objects = self.compiler.compile(sources,
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/ccompiler.py", line 574, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/home/sg/anaconda3/envs/shenfun/lib/python3.8/distutils/unixccompiler.py", line 120, in _compile
    raise CompileError(msg)
ImportError: Building module shenfun.optimization.cython.la failed: ["distutils.errors.CompileError: command '/home/sg/anaconda3/envs/shenfun/bin/x86_64-conda_cos6-linux-gnu-cc' failed with exit status 1\n"]

modifying dirichlet_poisson2D.py

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)

fix up docs in Getting started

I am looking at the docs and see that Function and Array also needs to be imported for the example to execute.

from shenfun import Basis, Function, Array
N = 8
T = Basis(N, 'Chebyshev', bc=None)

Also, later on in the code right before Operators you need to import Dx for this code to work.

From the Multipdimensional part, one needs to import TensorProductSpace, define comm and also import grad.

In coupled problems, the line beginning with SD does not work since int is not scriptable.

This is as far as I have gotten now. Will try and look at the rest another time.

Simple integration doing something funny

Hello Mikeal,

I am trying to learn how to integrate arrays to be able to compute energies and other such scalars.

Below is some simple code that I tried to integrate 1 to hopefully get one but I find that I actually get 0.712. Clearly this is not working how I thought. Could you clarify?

Also, I want to do this on a vector space and wonder if there is a way to use inner cleverly or simply write it out componet wise?

Thanks in advance,
Francis

    V1  = FunctionSpace(N[0], 'F', dtype='D', domain=(0, L[0]))
    tmp = Array(V1, val=1.0)
    l2_error = np.linalg.norm(tmp,ord=2)/N[0]

projecting with inhomogeneous boundary conditions

I am trying to compute the gradient of a function in a polynomial by Fourier space. It seems that if I sent the Dirichlet boundary conditions to be zero then I can do the projection easy enough, as is shown below. However, if I pick non-homogeneous boundary conditions I get an error.

Is there a way of doing this projection with inhomogeneous boundary conditions or should I instead do a decomposition into a line and a function that does have zero boundary conditions?

`#error

Traceback (most recent call last):
  File "qg_shenfun.py", line 86, in <module>
    ur = project(grad(pk),TV).backward(ur)
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/sitepackages/shenfun/forms/project.py", line 164, in project
    oa = b.solve(oa, oa)
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/sitepackages/shenfun/matrixbase.py", line 1235, in solve
    u = self.pmat.solve(b, u=u, axis=axis)
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/sitepackages/shenfun/matrixbase.py", line 625, in solve
    u = self.solver(b, u=u, axis=axis)
AttributeError: '_Legmatrix' object has no attribute 'solver'

# sample code

a = -1.
b =  1

# Function Spaces and grid
P0   = Basis(N[0], family='C', domain=(-L[0]/2,L[0]/2), bc=(a, b), scaled=True)
F1   = Basis(N[1], family='F', dtype='d', domain=(0., L[1]))
TF   = TensorProductSpace(comm, (P0, F1), axes=(0, 1))
TV   = VectorTensorProductSpace(TF)

X    = TF.local_mesh(True)

# Variables
pr, pk    = Array(TF), Function(TF)
ur, uk    = Array(TV), Function(TV)

pr[:] = np.sin(X[0])*np.sin(X[1])
pk  = TF.forward(pr, pk)
ur = project(grad(pk),TV).backward(ur)

KdV example has an error

I copied the KdV example into python and tried running it to find there was an error with the number of arguments in LinearRHS.

`# Error message

integrator.setup(dt)
  File "/home/fpoulin/software/anaconda3/envs/shenfun/lib/python3.8/site- 
packages/shenfun/utilities/integrators.py", line 195, in setup
    L = self.LinearRHS(**self.params)
 TypeError: LinearRHS() takes 0 positional arguments but 1 was given

`# KdV example

import numpy as np
from shenfun import *

N = 256
T = Basis(N, 'F', dtype='d')
u = TrialFunction(T)
v = TestFunction(T)
u_ = Array(T)
u_hat = Function(T)

def LinearRHS(**params):
    return -inner(Dx(u, 0, 3), v)

k = T.wavenumbers(scaled=True, eliminate_highest_freq=True)

def NonlinearRHS(u, u_hat, rhs, **params):
    rhs.fill(0)
    u_[:] = T.backward(u_hat, u_)
    rhs = T.forward(-0.5*u_**2, rhs)
    return rhs*1j*k   # return inner(grad(-0.5*Up**2), v)

A = 25.
B = 16.
x = T.points_and_weights()[0]
u_[:] = 3*A**2/np.cosh(0.5*A*(x-np.pi+2))**2 + 3*B**2/np.cosh(0.5*B*(x-np.pi+1))**2
u_hat = T.forward(u_, u_hat)

dt = 0.01/N**2
end_time = 0.006
integrator = ETDRK4(T, L=LinearRHS, N=NonlinearRHS)
integrator.setup(dt)
u_hat = integrator.solve(u_, u_hat, dt, (0, end_time))

pyFFTW runtimeerror with dirichlet_Poisson2D demo

With shenfun built and installed by conda running dirichlet_Poisson2D.py from demo folder produced

/home/miro3/miniconda2/envs/SHENFUN/lib/python2.7/site-packages/matplotlib/font_manager.py:280: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
'Matplotlib is building the font cache using fc-list. '

Traceback (most recent call last):
File "dirichlet_poisson2D.py", line 52, in
T = TensorProductSpace(comm, (SD, K1))
File "/home/miro3/miniconda2/envs/SHENFUN/lib/python2.7/site-packages/shenfun/tensorproductspace.py", line 79, in init
xfftn.plan(pencilB.subshape, axes, dtype, kw)
File "/home/miro3/miniconda2/envs/SHENFUN/lib/python2.7/site-packages/shenfun/chebyshev/bases.py", line 405, in plan
self.CT.plan(shape, axis, dtype, options)
File "/home/miro3/miniconda2/envs/SHENFUN/lib/python2.7/site-packages/shenfun/chebyshev/bases.py", line 147, in plan
xfftn_fwd = plan_fwd(U, axis=axis, **opts)
File "build/bdist.linux-x86_64/egg/pyfftw/builders/builders.py", line 549, in dct
File "build/bdist.linux-x86_64/egg/pyfftw/builders/_utils.py", line 205, in _Xfftn
File "pyfftw/pyfftw.pyx", line 1252, in pyfftw.pyfftw.FFTW.cinit (/home/mikael/anaconda2/conda-bld/work/pyfftw/pyfftw.c:11613)
RuntimeError: The data has an uncaught error that led to the planner returning NULL. This is a bug.

TypeError when applying inhomogeneous bcs in both x- and y- direction

Hi Mikael,
I am using shenfun to solve a vector Poisson problem on a 2D domain with inhomogeneous Dirichlet boundary conditions. Implementing inhomogeneous bcs in either x-direction or y-direction works perfectly fine. However, when I try to apply inhomogeneous bcs in both x-direction and y-direction I get a TypeError when I create the corresponding TensorProductSpace. I hope the code below shows the problem.

from mpi4py import MPI
from shenfun import FunctionSpace, TensorProductSpace
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD

family = 'legendre'
x, y = symbols("x,y")

# size of discretization
n = 20
N = [n, n+1]

# bases
BX_hom = x, y = symbols("x,y")(N[0], family=family, bc=(0, 0))
BY_hom = FunctionSpace(N[1], family=family, bc=(0, 0))
BX_inhom = FunctionSpace(N[0], family=family, bc=((1-y)**2*(1+y)**2, 0)) # inhomogeneous bc for v_x at x = 1
BY_inhom = FunctionSpace(N[1], family=family, bc=(0, (1-x)**2*(1+x)**2)) # inhomogeneous bc for v_x at y = 1

# spaces
subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T_hom = TensorProductSpace(subcomms, (BX_hom, BY_hom))

# working
T_inhom = TensorProductSpace(subcomms, (BX_hom, BY_inhom))

# non-working
# T_inhom = TensorProductSpace(subcomms, (BX_inhom, BY_inhom))

Best regards,
Thilo

spectraldns / shenfun Goto Github PK

shenfun's Introduction

Shenfun

Description

Installation

Dependencies

Contact

shenfun's People

Contributors

Stargazers

Watchers

Forkers

shenfun's Issues

Recommend Projects

Recommend Topics

Recommend Org