Cleanup tests (#254)

Clean up many tests.

dedalus/core/basis.py

@@ -550,8 +550,8 @@ def __mul__(self, other):
             if self.grid_params == other.grid_params:
                 # Everything matches except size, a, b
                 size = max(self.size, other.size)
-                a = max(self.a, other.a)
-                b = max(self.b, other.b)
+                a = self.a0
+                b = self.b0
                 return self.clone_with(size=size, a=a, b=b)
         if isinstance(other, SphereBasis):
             return other.__mul__(self)
@@ -566,7 +566,18 @@ def __matmul__(self, other):
 
     def __rmatmul__(self, other):
         # NCC (other) * operand (self)
-        return self.__mul__(other)
+        if other is None or other is self:
+            return self
+        if isinstance(other, Jacobi):
+            if self.grid_params == other.grid_params:
+                # Everything matches except size, a, b
+                size = max(self.size, other.size)
+                a = self.a
+                b = self.b
+                return self.clone_with(size=size, a=a, b=b)
+        if isinstance(other, SphereBasis):
+            return other.__mul__(self)
+        return NotImplemented
 
     # def include_mode(self, mode):
     #     return (0 <= mode < self.space.coeff_size)

dedalus/core/future.py

@@ -259,6 +259,11 @@ class FutureField(Future):
     """Class for deferred operations producing a Field."""
     future_type = Field
 
+    def __getitem__(self, layout):
+        """Evaluate and return data viewed in specified layout."""
+        field = self.evaluate()
+        return field[layout]
+
     @staticmethod
     def parse(string, namespace, dist):
         """Build FutureField from a string expression."""

dedalus/core/problems.py

@@ -44,6 +44,9 @@ def __init__(self, variables, namespace=None):
         # Build namespace via chainmap to react to upstream changes
         # Priority: local namespace, external namespace, parseables
         self.local_namespace = {}
+        for var in variables:
+            if var.name:
+                self.local_namespace[var.name] = var
         if namespace is None:
             self.namespace = ChainMap(self.local_namespace, parseables)
         else:
@@ -224,7 +227,8 @@ def __init__(self, *args, **kw):
         self.perturbations = [var.copy() for var in self.variables]
         for pert, var in zip(self.perturbations, self.variables):
             pert['c'] = 0
-            pert.name = 'δ'+var.name
+            if var.name:
+                pert.name = 'δ'+var.name
         self.LHS_variables = self.perturbations
 
     def _check_equation_conditions(self, eqn):

dedalus/core/solvers.py

@@ -8,6 +8,7 @@
 import scipy.linalg
 
 from . import subsystems
+from . import timesteppers
 from .evaluator import Evaluator
 from ..libraries.matsolvers import matsolvers
 from ..tools.config import config
@@ -483,7 +484,7 @@ class InitialValueSolver(SolverBase):
     ----------
     problem : Problem object
         Dedalus problem.
-    timestepper : Timestepper class
+    timestepper : Timestepper class or str
         Timestepper to use in evolving initial conditions.
     enforce_real_cadence : int, optional
         Iteration cadence for enforcing Hermitian symmetry on real variables (default: 100).
@@ -535,6 +536,8 @@ def __init__(self, problem, timestepper, enforce_real_cadence=100, warmup_iterat
         F_handler.build_system()
         self.F = F_handler.fields
         # Initialize timestepper
+        if isinstance(timestepper, str):
+            timestepper = timesteppers.schemes[timestepper]
         self.timestepper = timestepper(self)
         # Attributes
         self.sim_time = self.initial_sim_time = problem.time.allreduce_data_max(layout='g')

dedalus/extras/flow_tools.py

@@ -198,7 +198,7 @@ def compute_timestep(self):
             if (iteration-1) <= self.solver.initial_iteration:
                 return self.stored_dt
             # Sum across frequencies for each local grid point
-            local_freqs = np.sum(np.abs(field['g']) for field in self.frequencies.fields.values())
+            local_freqs = sum(np.abs(field['g']) for field in self.frequencies.fields.values())
             # Compute new timestep from max frequency across all grid points
             max_global_freq = self.reducer.global_max(local_freqs)
             if max_global_freq == 0.:

dedalus/tests/ball_diffusion_analytical_eigenvalues.py

@@ -1,11 +1,12 @@
+"""Laplacian eigenvalues in the ball with various boundary conditions."""
+
 import numpy as np
 from scipy.special import spherical_jn as j
 
-def dispersion_zeros(ell,n,a=0,r0=1,guess=None,imax=20,nk=10,eps=0.1):
 
+def dispersion_zeros(ell, n, a=0, r0=1, guess=None, imax=20, nk=10, eps=0.1):
     def F(k,deriv=False):
-        return j(ell,k*r0,derivative=deriv) - a*j(ell+2,k*r0,derivative=deriv)
-
+        return j(ell, k*r0, derivative=deriv) - a*j(ell+2, k*r0, derivative=deriv)
     if guess == None:
         kmax = np.pi*(n+ell/2 + eps)/r0
         k = np.linspace(0,kmax,int(kmax*nk))
@@ -14,46 +15,42 @@ def F(k,deriv=False):
         k = 0.5*(k[i]+k[i+1])
     else:
         k = guess
-
     for i in range(imax):
         dk =  F(k)/F(k,deriv=True)
         k -= dk
-
     print('dk =',np.max(np.abs(dk)))
-
     return k
 
 
-def wavenumbers(ell,n,BC, **kwargs):
-
-    k = {'toroidal':0,'poloidal':0}
-
+def wavenumbers(ell, n, BC, **kw):
+    k = {'toroidal':0, 'poloidal':0}
     if BC=='Bessel':
-        k = dispersion_zeros(ell,n,**kwargs)
+        k = dispersion_zeros(ell, n, **kw)
     elif BC=="no-slip":
-        k['toroidal'] = dispersion_zeros(ell,n,**kwargs)
-        k['poloidal'] = dispersion_zeros(ell+1,n,**kwargs)
+        k['toroidal'] = dispersion_zeros(ell, n, **kw)
+        k['poloidal'] = dispersion_zeros(ell+1, n, **kw)
     elif BC=="stress-free":
         if ell == 1:
-            k['toroidal'] = dispersion_zeros(2,n,**kwargs)
+            k['toroidal'] = dispersion_zeros(2, n, **kw)
         else:
-            k['toroidal'] = dispersion_zeros(ell-1,n,a=(ell+2)/(ell-1),**kwargs)
-        k['poloidal'] = dispersion_zeros(ell,n,a=2/(2*ell+1),**kwargs)
+            k['toroidal'] = dispersion_zeros(ell-1, n, a=(ell+2)/(ell-1), **kw)
+        k['poloidal'] = dispersion_zeros(ell, n, a=2/(2*ell+1), **kw)
     elif BC=="potential":
-        k['toroidal'] = dispersion_zeros(ell-1,n,**kwargs)
-        k['poloidal'] = dispersion_zeros(ell,n,**kwargs)
+        k['toroidal'] = dispersion_zeros(ell-1, n, **kw)
+        k['poloidal'] = dispersion_zeros(ell, n, **kw)
     elif BC=="conducting":
-        k['toroidal'] = dispersion_zeros(ell,n,**kwargs)
-        k['poloidal'] = dispersion_zeros(ell-1,n,a=ell/(ell+1),**kwargs)
+        k['toroidal'] = dispersion_zeros(ell, n, **kw)
+        k['poloidal'] = dispersion_zeros(ell-1, n, a=ell/(ell+1), **kw)
     elif BC=="pseudo":
-        k['toroidal'] = dispersion_zeros(ell-1,n,a=ell/(ell+1),**kwargs)
-        k['poloidal'] = dispersion_zeros(ell,n,**kwargs)
+        k['toroidal'] = dispersion_zeros(ell-1, n, a=ell/(ell+1), **kw)
+        k['poloidal'] = dispersion_zeros(ell, n, **kw)
     else:
         raise ValueError("BC='{}' is not a valid choice".format(BC))
-
     return k
 
-def eigenvalues(ell,n,BC, **kwargs):
-    k = wavenumbers(ell,n,BC, **kwargs)
-    k = np.sort(np.concatenate((k['toroidal'],k['poloidal'])))
+
+def eigenvalues(ell, n, BC, **kw):
+    k = wavenumbers(ell, n, BC, **kw)
+    k = np.sort(np.concatenate((k['toroidal'], k['poloidal'])))
     return k**2
+

dedalus/tests/test_cartesian_ncc.py

@@ -1,4 +1,5 @@
-"""Jacobi and Fourier NCC tests."""
+"""Test Fourier and Jacobi NCCs."""
+# TODO: add tests inovlving constant NCCs/operands
 
 import pytest
 import numpy as np
@@ -27,10 +28,9 @@ def build_jacobi(N, a0, b0, bounds, dealias, dtype):
 
 
 @pytest.mark.parametrize('N', [16])
-@pytest.mark.parametrize('a0', [-1/2])
-@pytest.mark.parametrize('b0', [-1/2])
-@pytest.mark.parametrize('k_ncc', [0, 1, 2])
-@pytest.mark.parametrize('k_arg', [0, 1, 2])
+@pytest.mark.parametrize('a0, b0', [(-1/2, -1/2), (0, 0)])
+@pytest.mark.parametrize('k_ncc', [0, 1])
+@pytest.mark.parametrize('k_arg', [0, 1])
 @pytest.mark.parametrize('dealias', [3/2])
 @pytest.mark.parametrize('dtype', [np.float64, np.complex128])
 def test_eval_jacobi_ncc(N, a0, b0, k_ncc, k_arg, dealias, dtype):
@@ -51,12 +51,9 @@ def test_eval_jacobi_ncc(N, a0, b0, k_ncc, k_arg, dealias, dtype):
     w1.store_ncc_matrices(vars, solver.subproblems)
     w0 = w0.evaluate()
     w1 = w1.evaluate_as_ncc()
-    # Remove last kmax coeffs which have dealias-before-conversion errors
-    k_max = max(k_ncc, k_arg)
-    if k_max > 0:
-        w0['c'][-2*k_max:] = 0
-        w1['c'][-2*k_max:] = 0
-    assert np.allclose(w0['c'], w1['c'])
+    # Remove last 2*k_arg coeffs which have dealias-before-conversion errors
+    w2 = (w1 - w0).evaluate()
+    assert np.allclose(w2['c'][:-2*k_arg], 0)
 
 
 @pytest.mark.parametrize('N', [16])
@@ -84,8 +81,7 @@ def test_eval_fourier_ncc(N, dealias, dtype):
 
 
 @pytest.mark.parametrize('N', [16])
-@pytest.mark.parametrize('a0', [-1/2])
-@pytest.mark.parametrize('b0', [-1/2])
+@pytest.mark.parametrize('a0, b0', [(-1/2, -1/2)])
 @pytest.mark.parametrize('f_rank', [0, 1])
 @pytest.mark.parametrize('g_rank', [0, 1])
 @pytest.mark.parametrize('dealias', [3/2])
@@ -117,11 +113,10 @@ def test_eval_fourier_jacobi_ncc(N, a0, b0, f_rank, g_rank, dealias, dtype):
 
 
 @pytest.mark.parametrize('N', [16])
-@pytest.mark.parametrize('a0', [-1/2, 0])
-@pytest.mark.parametrize('b0', [-1/2, 0])
-@pytest.mark.parametrize('k_ncc', [0, xfail_param(1, "dealias before truncation", run=False)])
-@pytest.mark.parametrize('k_arg', [0, xfail_param(1, "dealias before truncation", run=False)])
-@pytest.mark.parametrize('dealias', [3/2, 2])
+@pytest.mark.parametrize('a0, b0', [(-1/2, -1/2), (0, 0)])
+@pytest.mark.parametrize('k_ncc', [0, 1])
+@pytest.mark.parametrize('k_arg', [0, xfail_param(1, "dealias before truncation", run=True)])
+@pytest.mark.parametrize('dealias', [3/2])
 @pytest.mark.parametrize('dtype', [np.float64, np.complex128])
 def test_solve_jacobi_ncc(N, a0, b0, k_ncc, k_arg, dealias, dtype):
     """Solve f(x) * u(x) = f(x) * g(x)."""
@@ -141,7 +136,7 @@ def test_solve_jacobi_ncc(N, a0, b0, k_ncc, k_arg, dealias, dtype):
 
 
 @pytest.mark.parametrize('N', [16])
-@pytest.mark.parametrize('dealias', [3/2, 2])
+@pytest.mark.parametrize('dealias', [3/2])
 @pytest.mark.parametrize('dtype', [np.float64, np.complex128])
 def test_solve_fourier_ncc(N, dealias, dtype):
     """Solve f(x) * u(x) = f(x) * g(x)."""
@@ -163,97 +158,3 @@ def test_solve_fourier_ncc(N, dealias, dtype):
     g.change_scales(1)
     assert np.allclose(u['g'], g['g'])
 
-
-# Lx=2
-# Ly=2
-# Lz=1
-# @CachedMethod
-# def build_3d_box(Nx, Ny, Nz, dealias, dtype, k=0):
-#     c = coords.CartesianCoordinates('x', 'y', 'z')
-#     d = distributor.Distributor((c,))
-#     if dtype == np.complex128:
-#         xb = basis.ComplexFourier(c.coords[0], size=Nx, bounds=(0, Lx))
-#         yb = basis.ComplexFourier(c.coords[1], size=Ny, bounds=(0, Ly))
-#     elif dtype == np.float64:
-#         xb = basis.RealFourier(c.coords[0], size=Nx, bounds=(0, Lx))
-#         yb = basis.RealFourier(c.coords[1], size=Ny, bounds=(0, Ly))
-#     zb = basis.ChebyshevT(c.coords[2], size=Nz, bounds=(0, Lz))
-#     b = (xb, yb, zb)
-#     x = xb.local_grid(1)
-#     y = yb.local_grid(1)
-#     z = zb.local_grid(1)
-#     return c, d, b, x, y, z
-
-# @CachedMethod
-# def build_2d_box(Nx, Nz, dealias, dtype, k=0):
-#     c = coords.CartesianCoordinates('x', 'z')
-#     d = distributor.Distributor((c,))
-#     if dtype == np.complex128:
-#         xb = basis.ComplexFourier(c.coords[0], size=Nx, bounds=(0, Lx))
-#     elif dtype == np.float64:
-#         xb = basis.RealFourier(c.coords[0], size=Nx, bounds=(0, Lx))
-#     zb = basis.ChebyshevT(c.coords[1], size=Nz, bounds=(0, Lz))
-#     b = (xb, zb)
-#     x = xb.local_grid(1)
-#     z = zb.local_grid(1)
-#     return c, d, b, x, z
-
-# Nx_range = [8]
-# Ny_range = [8]
-# Nz_range = [8]
-# dealias_range = [1, 3/2]
-
-# @pytest.mark.parametrize('Nx', Nx_range)
-# @pytest.mark.parametrize('Ny', Ny_range)
-# @pytest.mark.parametrize('Nz', Nz_range)
-# @pytest.mark.parametrize('dealias', dealias_range)
-# @pytest.mark.parametrize('basis', [build_3d_box])
-# @pytest.mark.parametrize('dtype', [np.complex128, np.float64])
-# def test_3d_scalar_ncc_scalar(Nx, Ny, Nz, dealias, basis, dtype):
-#     c, d, b, x, y, z = basis(Nx, Ny, Nz, dealias, dtype)
-#     u = field.Field(dist=d, bases=b, dtype=dtype)
-#     ncc = field.Field(name='ncc', dist=d, bases=(b[-1],), dtype=dtype)
-#     ncc['g'] = 1/(1+z**2)
-#     problem = problems.LBVP([u])
-#     problem.add_equation((ncc*u, 1))
-#     # Solver
-#     solver = solvers.LinearBoundaryValueSolver(problem)
-#     solver.solve()
-#     assert np.allclose(u['g'], 1/ncc['g'])
-
-# @pytest.mark.parametrize('Nx', Nx_range)
-# @pytest.mark.parametrize('Ny', Ny_range)
-# @pytest.mark.parametrize('Nz', Nz_range)
-# @pytest.mark.parametrize('dealias', dealias_range)
-# @pytest.mark.parametrize('basis', [build_3d_box])
-# @pytest.mark.parametrize('dtype', [np.complex128, np.float64])
-# def test_3d_scalar_ncc_scalar_separable(Nx, Ny, Nz, dealias, basis, dtype):
-#     c, d, b, x, y, z = basis(Nx, Ny, Nz, dealias, dtype)
-#     u = field.Field(dist=d, bases=(b[0], b[1]), dtype=dtype)
-#     v = field.Field(dist=d, bases=b, dtype=dtype)
-#     ncc = field.Field(name='ncc', dist=d, bases=(b[-1],), dtype=dtype)
-#     ncc['g'] = 1/(1+z**2)
-#     problem = problems.LBVP([v, u])
-#     problem.add_equation((v, 0)) # Hack by adding v since problem must have full dimension
-#     problem.add_equation(((ncc*u)(z=1), ncc(z=1)))
-#     # Solver
-#     solver = solvers.LinearBoundaryValueSolver(problem)
-#     solver.solve()
-#     assert np.allclose(u['g'], 1)
-
-# @pytest.mark.parametrize('Nx', Nx_range)
-# @pytest.mark.parametrize('Nz', Nz_range)
-# @pytest.mark.parametrize('dealias', dealias_range)
-# @pytest.mark.parametrize('basis', [build_2d_box])
-# @pytest.mark.parametrize('dtype', [np.complex128, np.float64])
-# def test_implicit_transpose_2d_tensor(Nx, Nz, dealias, basis, dtype):
-#     c, d, b, x, z = basis(Nx, Nz, dealias, dtype)
-#     u = field.Field(dist=d, bases=b, dtype=dtype)
-#     ncc = field.Field(name='ncc', dist=d, bases=(b[-1],), dtype=dtype)
-#     ncc['g'] = 1/(1+z**2)
-#     problem = problems.LBVP([u])
-#     problem.add_equation((ncc*u, 1))
-#     # Solver
-#     solver = solvers.LinearBoundaryValueSolver(problem)
-#     solver.solve()
-#     assert np.allclose(u['g'], 1/ncc['g'])