Some minor fixes and bump version

demo/ChannelFlow2D.py

@@ -90,10 +90,10 @@ def __init__(self,
         self.TD = TensorProductSpace(comm, (self.D0, self.F1), modify_spaces_inplace=True) # Streamwise velocity
         self.TC = TensorProductSpace(comm, (self.C0, self.F1), modify_spaces_inplace=True) # No bc
         self.BD = VectorSpace([self.TB, self.TD])  # Velocity vector space
-        self.CD = VectorSpace(self.TD)             # Dirichlet vector space
+        self.CD = VectorSpace(self.TC)             # Dirichlet vector space
 
         # Padded space for dealiasing
-        self.TDp = self.TD.get_dealiased(padding_factor)
+        self.TDp = self.TC.get_dealiased(padding_factor)
 
         self.u_ = Function(self.BD)      # Velocity vector solution
         self.H_ = Function(self.CD)      # convection
@@ -110,9 +110,9 @@ def __init__(self,
         # Classes for fast projections. All are not used except if self.conv=0
         self.dudx = Project(Dx(self.u_[0], 0, 1), self.TD)
         if self.conv == 0:
-            self.dudy = Project(Dx(self.u_[0], 1, 1), self.TB)
+            self.dudy = Project(Dx(self.u_[0], 1, 1), self.TC)
             self.dvdx = Project(Dx(self.u_[1], 0, 1), self.TC)
-            self.dvdy = Project(Dx(self.u_[1], 1, 1), self.TD)
+            self.dvdy = Project(Dx(self.u_[1], 1, 1), self.TC)
 
         self.curl = Project(curl(self.u_), self.TC)
         self.divu = Project(div(self.u_), self.TD)
@@ -130,8 +130,8 @@ def __init__(self,
         v = TestFunction(self.TB)
 
         # Chebyshev matrices are not sparse, so need a tailored solver. Legendre has simply 5 nonzero diagonals and can use generic solvers.
-        #sol1 = chebyshev.la.Biharmonic if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
-        sol1 = la.SolverGeneric1ND
+        sol1 = chebyshev.la.Biharmonic if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
+        #sol1 = la.SolverGeneric1ND
 
         self.pdes = {
 
@@ -151,7 +151,7 @@ def __init__(self,
         # v. Momentum equation for Fourier wavenumber 0
         if comm.Get_rank() == 0:
             v0 = TestFunction(self.D00)
-            self.h1 = Function(self.D00)  # Copy from H_[1, :, 0, 0] (cannot use view since not contiguous)
+            self.h1 = Function(self.C00)  # Copy from H_[1, :, 0, 0] (cannot use view since not contiguous)
             source = Array(self.C00)
             source[:] = -self.dpdy        # dpdy set by subclass
             sol = chebyshev.la.Helmholtz if self.B0.family() == 'chebyshev' else la.Solver
@@ -164,13 +164,14 @@ def __init__(self,
                           solver=sol,
                           latex=r"\frac{\partial v}{\partial t} = \nu \frac{\partial^2 v}{\partial x^2} - N_y - \frac{\partial p}{\partial y}"),
             }
-
-    def convection(self):
+    
+    def convection(self, rk):
         H = self.H_.v
         self.up = self.u_.backward(padding_factor=self.padding_factor)
         up = self.up.v
         if self.conv == 0:
-            dudxp = self.dudx().backward(padding_factor=self.padding_factor).v
+            dudx = self.dudx() if rk == 0 else self.dudx.output_array
+            dudxp = dudx.backward(padding_factor=self.padding_factor).v
             dudyp = self.dudy().backward(padding_factor=self.padding_factor).v
             dvdxp = self.dvdx().backward(padding_factor=self.padding_factor).v
             dvdyp = self.dvdy().backward(padding_factor=self.padding_factor).v
@@ -192,7 +193,7 @@ def compute_v(self, rk):
 
         # Find velocity components v from div. constraint
         u[1] = 1j*self.dudx()/self.K[1]
-
+            
         # Still have to compute for wavenumber = 0, 0
         if comm.Get_rank() == 0:
             # v component
@@ -206,7 +207,6 @@ def compute_pressure(self):
             self.d2udx2 = Project(self.nu*Dx(self.u_[0], 0, 2), self.TC)
             N0 = self.N0 = FunctionSpace(self.N[0], self.B0.family(), bc={'left': {'N': self.d2udx2()}, 'right': {'N': self.d2udx2()}})
             TN = self.TN = TensorProductSpace(comm, (N0, self.F1), modify_spaces_inplace=True)
-            #sol = chebyshev.la.Helmholtz if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
             sol = la.SolverGeneric1ND
             self.divH = Inner(TestFunction(TN), -div(self.H_))
             self.solP = sol(inner(TestFunction(TN), div(grad(TrialFunction(TN)))))
@@ -251,15 +251,15 @@ def tofile(self, tstep):
         self.file_u.write(tstep, {'u': [self.u_.backward(mesh='uniform')]}, as_scalar=True)
 
     def prepare_step(self, rk):
-        self.convection()
+        self.convection(rk)
 
     def assemble(self):
         for pde in self.pdes.values():
             pde.assemble()
         if comm.Get_rank() == 0:
             for pde in self.pdes1d.values():
                 pde.assemble()
-
+    
     def solve(self, t=0, tstep=0, end_time=1000):
         self.assemble()
         while t < end_time-1e-8:

demo/KleinGordon.py

@@ -29,7 +29,7 @@
 ue = 0.1*exp(-(x**2 + y**2 + z**2))
 
 # Size of discretization
-N = (32, 32, 32)
+N = (31, 32, 33)
 
 # Defocusing or focusing
 gamma = 1
@@ -104,8 +104,8 @@ def update(self, fu, fu_hat, t, tstep, **params):
     if rank == 0 and tstep % params['plot_tstep'] == 0 and params['plot_tstep'] > 0:
         fu = fu_hat.backward(fu)
         f, u = fu[:]
-        image.ax.clear()
-        image.ax.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100)
+        #image.axes.clear()
+        image.axes.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100)
         plt.pause(1e-6)
         transformed = True
 
@@ -144,7 +144,7 @@ def update(self, fu, fu_hat, t, tstep, **params):
            'write_tstep': (50, {'fu': [fu]}),
            'Compute_energy': 100,
            'plot_tstep': 100,
-           'end_time': 100.,
+           'end_time': 10.,
            'file': file0}
     dt = 0.005
     #integrator = ETDRK4(TT, N=NonlinearRHS, update=update, **par)

demo/OrrSommerfeld.py

@@ -82,10 +82,8 @@ def plot(self, t, tstep):
             if comm.Get_rank() == 0 and comm.Get_size() == 1:
                 ub = self.u_.backward(self.ub)
                 X = self.X
-                self.im1.axes.clear()
                 self.im1.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[0, :, :, 0], 100)
                 self.im1.autoscale()
-                self.im2.axes.clear()
                 self.im2.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[1, :, :, 0], 100)
                 self.im2.autoscale()
                 self.im3.set_UVC(ub[1, :, :, 0]-(1-self.X[0][:, :, 0]**2), ub[0, :, :, 0])

demo/OrrSommerfeld2D.py

@@ -5,27 +5,26 @@
 
 class OrrSommerfeld(KMM):
 
-    def __init__(self, N=(32, 32), domain=((-1, 1), (0, 2*np.pi)), Re=8000.,
+    def __init__(self, N=(32, 32), domain=((-1, 1), (0, 2*np.pi)), Re=8000., alfa=1.,
                  dt=0.1, conv=0, modplot=100, modsave=1e8, moderror=100, filename='KMM',
                  family='C', padding_factor=(1, 1.5), checkpoint=1000, timestepper='IMEXRK3'):
-        KMM.__init__(self, N=N, domain=domain, nu=1/Re, dt=dt, conv=conv, modplot=modplot,
+        KMM.__init__(self, N=N, domain=(domain[0], (domain[1][0], domain[1][1]/alfa)), nu=1/Re, dt=dt, conv=conv, modplot=modplot,
                      modsave=modsave, moderror=moderror, filename=filename, family=family,
                      padding_factor=padding_factor, checkpoint=checkpoint, timestepper=timestepper,
                      dpdy=-2/Re)
         self.Re = Re
+        self.alfa = alfa
 
     def initialize(self, from_checkpoint=False):
         if from_checkpoint:
             return self.init_from_checkpoint()
 
         from OrrSommerfeld_eigs import OrrSommerfeld
-        self.OS = OS = OrrSommerfeld(Re=self.Re, N=128)
+        self.OS = OS = OrrSommerfeld(Re=self.Re, N=128, alfa=self.alfa)
         eigvals, eigvectors = OS.solve(False)
         OS.eigvals, OS.eigvectors = eigvals, eigvectors
         self.initOS(OS, eigvals, eigvectors, self.ub)
         self.u_ = self.ub.forward(self.u_)
-        # Compute convection from data in context (i.e., context.U_hat and context.g)
-        # This is the convection at t=0
         self.e0 = 0.5*dx(self.ub[0]**2+(self.ub[1]-(1-self.X[0]**2))**2)
         self.acc = np.zeros(1)
         return 0, 0
@@ -37,16 +36,16 @@ def initOS(self, OS, eigvals, eigvectors, U, t=0.):
         OS.eigval = eigval
         for j in range(U.shape[2]):
             y = X[1][0, j]
-            v = (1-x**2) + 1e-7*np.real(dphidy*np.exp(1j*(y-eigval*t)))
-            u = -1e-7*np.real(1j*phi*np.exp(1j*(y-eigval*t)))
+            v = (1-x**2) + 1e-7*np.real(dphidy*np.exp(1j*self.alfa*(y-eigval*t)))
+            u = -1e-7*np.real(1j*self.alfa*phi*np.exp(1j*self.alfa*(y-eigval*t)))
             U[0, :, j] = u
             U[1, :, j] = v
 
     def compute_error(self, t):
         ub = self.u_.backward(self.ub)
         pert = (ub[1] - (1-self.X[0]**2))**2 + ub[0]**2
         e1 = 0.5*dx(pert)
-        exact = np.exp(2*np.imag(self.OS.eigval)*t)
+        exact = np.exp(2*np.imag(self.alfa*self.OS.eigval)*t)
         U0 = self.work[(ub, 0, True)]
         self.initOS(self.OS, self.OS.eigvals, self.OS.eigvectors, U0, t=t)
         pert = (ub[0] - U0[0])**2 + (ub[1] - U0[1])**2
@@ -80,10 +79,8 @@ def plot(self, t, tstep):
             if comm.Get_rank() == 0 and comm.Get_size() == 1:
                 ub = self.u_.backward(self.ub)
                 X = self.X
-                self.im1.axes.clear()
                 self.im1.axes.contourf(X[1], X[0], ub[0], 100)
                 self.im1.autoscale()
-                self.im2.axes.clear()
                 self.im2.axes.contourf(X[1], X[0], ub[1], 100)
                 self.im2.autoscale()
                 self.im3.set_UVC(ub[1]-(1-self.X[0]**2), ub[0])
@@ -105,25 +102,26 @@ def print_energy_and_divergence(self, t, tstep):
     from mpi4py_fft import generate_xdmf
     t0 = time()
     N = (128, 32)
-    config['optimization']['mode'] = 'numba'
+    config['optimization']['mode'] = 'cython'
     d = {
         'N': N,
         'Re': 8000.,
-        'dt': 0.01,
+        'dt': 0.002,
         'filename': f'KMM_OS_{N[0]}_{N[1]}',
         'conv': 0,
-        'modplot': 100,
+        'alfa': 1,
+        'modplot': -1,
         'modsave': 1000,
         'moderror': 10,
         'family': 'C',
         'checkpoint': 10000000,
         'padding_factor': 1,
-        'timestepper': 'IMEXRK222'
+        'timestepper': 'IMEXRK443'
         }
     OS = True
     c = OrrSommerfeld(**d)
     t, tstep = c.initialize(from_checkpoint=False)
-    c.solve(t=t, tstep=tstep, end_time=1)
+    c.solve(t=t, tstep=tstep, end_time=0.1)
     print('Computing time %2.4f'%(time()-t0))
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')

shenfun/__init__.py

@@ -29,7 +29,7 @@
 """
 #pylint: disable=wildcard-import,no-name-in-module
 
-__version__ = '4.2.0'
+__version__ = '4.2.1'
 __author__ = 'Mikael Mortensen'
 
 import numpy as np

shenfun/forms/arguments.py

@@ -346,7 +346,7 @@ def FunctionSpace(N, family='Fourier', bc=None, dtype='d', quad=None,
             else:
                 raise NotImplementedError
             B = bases[''.join(key)]
-
+        
         return B(N, **par)
 
     elif family.lower() in ('hermite', 'h'):

shenfun/la.py

@@ -132,6 +132,10 @@ def perform_lu(self):
         if self._lu is None:
             if isinstance(self.mat, SparseMatrix):
                 self.mat = self.mat.diags('csc')
+            self.mat.eliminate_zeros()
+            if self.mat.nnz < 2:
+                self._lu = None
+                return self._lu
             self._lu = splu(self.mat, permc_spec=config['matrix']['sparse']['permc_spec'])
             self.dtype = self.mat.dtype.char
             self._inner_arg = (self._lu, self.dtype)
@@ -238,7 +242,8 @@ def __call__(self, b, u=None, axis=0, constraints=()):
         b = self.apply_bcs(b, u, axis=axis)
         b = self.apply_constraints(b, constraints, axis=axis)
         lu = self.perform_lu() # LU must be performed after constraints, because constraints modify the matrix
-        u = self.solve(b, u, axis, lu)
+        if not lu is None:
+            u = self.solve(b, u, axis, lu)
         if hasattr(u, 'set_boundary_dofs'):
             u.set_boundary_dofs()
         return u
@@ -838,7 +843,6 @@ def __init__(self, mats, format=None):
         SparseMatrixSolver.__init__(self, mats)
         self.mat = self.mat.diags(format)
 
-
 class SolverGeneric2ND:
     """Generic solver for problems consisting of tensorproduct matrices
     containing two non-diagonal submatrices.
@@ -1371,15 +1375,17 @@ def solve(self, u, b, solvers1D, naxes):
             for i, sol in enumerate(solvers1D):
                 s[paxes[0]] = i
                 s0 = tuple(s)
-                sol.inner_solve(u[s0], sol._inner_arg)
+                if sol._inner_arg is not None:
+                    sol.inner_solve(u[s0], sol._inner_arg)
 
         elif u.ndim == 3:
             for i, m in enumerate(solvers1D):
                 s[paxes[0]] = i
                 for j, sol in enumerate(m):
                     s[paxes[1]] = j
                     s0 = tuple(s)
-                    sol.inner_solve(u[s0], sol._inner_arg)
+                    if sol._inner_arg is not None:
+                        sol.inner_solve(u[s0], sol._inner_arg)
 
     def __call__(self, b, u=None, constraints=(), fast=True):
         """Solve problem with one non-diagonal direction
@@ -1472,9 +1478,7 @@ def __call__(self, b, u=None, constraints=()):
         space = b.function_space()
         if u is None:
             u = Function(space)
-        else:
-            assert u.shape == b.shape
-
+
         if self.bc_mat: # Add contribution to right hand side due to inhomogeneous boundary conditions
             u.set_boundary_dofs()
             w0 = np.zeros_like(b)

shenfun/matrixbase.py

@@ -446,6 +446,9 @@ def get_solver(self):
         """
         from .la import (Solve, TDMA, TDMA_O, FDMA, TwoDMA, ThreeDMA, PDMA,
             DiagMA, HeptaDMA)
+        if self.scale == 0:
+            return Solve
+
         if len(self) == 1:
             if list(self.keys())[0] == 0:
                 return DiagMA

shenfun/utilities/findbasis.py

@@ -46,9 +46,10 @@ def get_stencil_matrix(bcs, family, alpha=None, beta=None, gn=1):
     r = []
     for key in bcs.orderednames():
         k, v = key[0], key[1:]
-        if v == 'R':
-            alfa = bcs[lra[k]]['R'][0]
-            f = [bnd_values(k=0)[lr[k]], bnd_values(k=1)[lr[k]]]
+        if v in 'WR':
+            k0 = 0 if v == 'R' else 1
+            alfa = bcs[lra[k]][v][0]
+            f = [bnd_values(k=k0)[lr[k]], bnd_values(k=k0+1)[lr[k]]]
             s.append([sp.simplify(f[0](n+j)+alfa*f[1](n+j)) for j in range(1, 1+bcs.num_bcs())])
             r.append(-sp.simplify(f[0](n)+alfa*f[1](n)))
         else:
@@ -100,9 +101,10 @@ def _computematrix(first):
         s = []
         for key in bcs.orderednames():
             k, v = key[0], key[1:]
-            if v == 'R':
-                alfa = bcs[lra[k]]['R'][0]
-                f = [bnd_values(k=0)[lr[k]], bnd_values(k=1)[lr[k]]]
+            if v in 'WR':
+                k0 = 0 if v == 'R' else 1
+                alfa = bcs[lra[k]][v][0]
+                f = [bnd_values(k=k0)[lr[k]], bnd_values(k=k0+1)[lr[k]]]
                 s.append([sp.simplify(f[0](j)+alfa*f[1](j)) for j in range(first, first+bcs.num_bcs())])
             else:
                 f = bnd_values(k=bc[v])[lr[k]]