Function test_ball_bessel_eigenfunction

dedalus/tests/test_evp.py:153–193

Test eigenfunctions of Bessel equation in ball bases.

(Nphi, Ntheta, Nr, radius, alpha, dtype, ell)

@pytest.mark.parametrize('dtype', [np.complex128])
@pytest.mark.parametrize('ell', [0, 1, 2, 3])
def test_ball_bessel_eigenfunction(Nphi, Ntheta, Nr, radius, alpha, dtype, ell):
    """Test eigenfunctions of Bessel equation in ball bases."""
    # Bases
    c = d3.SphericalCoordinates('phi', 'theta', 'r')
    d = d3.Distributor(c, dtype=dtype)
    b = d3.BallBasis(c, (Nphi, Ntheta, Nr), radius=radius, alpha=alpha, dtype=dtype)
    phi, theta, r = d.local_grids(b, scales=(1, 1, 1))
    # Fields
    u = d.Field(bases=b)
    tau = d.Field(bases=b.surface)
    s = d.Field()
    # Problem
    lift = lambda A: d3.Lift(A, b, -1)
    problem = d3.EVP([u, tau], s, namespace=locals())
    problem.add_equation("s*u + lap(u) + lift(tau) = 0")
    problem.add_equation("u(r=radius) = 0")
    # Solver
    solver = problem.build_solver()
    # TODO: clean this up with group selection interface
    for sp in solver.subproblems:
        if sp.group[1] == ell:
            break
    else:
        raise ValueError(f"Could not find subproblem with ell = {ell}")
    solver.solve_dense(sp)
    # Compare eigenfunction
    i_sort = np.argsort(solver.eigenvalues)
    solver.eigenvalues = solver.eigenvalues[i_sort]
    solver.eigenvectors = solver.eigenvectors[:,i_sort]
    for Neig in [0, 1, 5, 10]:
        solver.set_state(Neig, sp.subsystems[0]) # m = 0 mode
        u.change_layout(d.layouts[1])
        local_m, local_ell, local_n = u.layout.local_group_arrays(u.domain, u.scales)
        radial_eigenfunction = u.data[(local_m == 0)*(local_ell == ell)]
        i_max = np.argmax(np.abs(radial_eigenfunction))
        radial_eigenfunction /= radial_eigenfunction[i_max]
        k = np.sqrt(solver.eigenvalues[Neig])
        sol = spec.jv(ell+1/2, k*r) / np.sqrt(k*r)
        sol = sol.ravel()
        sol /= sol[i_max]
        assert np.allclose(radial_eigenfunction, sol)


@pytest.mark.parametrize('Nphi', [8])