Function test_ball_diffusion

dedalus/tests/test_evp.py:204–252

Test eigenvalues of vector diffusion equation in ball bases with various BCs.

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

@pytest.mark.parametrize('ell', [1, 2, 3])
@pytest.mark.parametrize('bc', ['no-slip', 'stress-free', 'potential', 'conducting', 'pseudo'])
def test_ball_diffusion(Nphi, Ntheta, Nr, radius, alpha, dtype, ell, bc):
    """Test eigenvalues of vector diffusion equation in ball bases with various BCs."""
    # 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)
    # Fields
    φ = d.Field(bases=b)
    A = d.VectorField(c, bases=b)
    tA = d.VectorField(c, bases=b.surface)
    λ = d.Field()
    # Problem
    lift = lambda A: d3.Lift(A, b, -1)
    problem = d3.EVP([φ, A, tA], λ, namespace=locals())
    problem.add_equation("div(A) = 0")
    problem.add_equation("-λ*A + grad(φ) - lap(A) + lift(tA) = 0")
    if bc == 'no-slip':
        problem.add_equation("A(r=radius) = 0")
    elif bc == 'stress-free':
        problem.namespace['E'] = (d3.grad(A) + d3.grad(A).T) / 2
        problem.add_equation("radial(A(r=radius)) = 0")
        problem.add_equation("angular(radial(E(r=radius))) = 0")
    elif bc == 'potential':
        ell_func = lambda ell: ell+1
        problem.namespace['ell_1'] = lambda A: d3.SphericalEllProduct(A, c, ell_func)
        problem.add_equation("radial(grad(A)(r=radius)) + ell_1(A)(r=radius)/radius = 0")
    elif bc == 'conducting':
        problem.add_equation("φ(r=radius) = 0")
        problem.add_equation("angular(A(r=radius)) = 0")
    elif bc == 'pseudo':
        problem.add_equation("radial(A(r=radius)) = 0")
        problem.add_equation("angular(curl(A)(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 eigenvalues
    i_sort = np.argsort(solver.eigenvalues)
    solver.eigenvalues = solver.eigenvalues[i_sort]
    λ_analytic = analytic_eigenvalues(ell, Nr, bc, r0=radius)
    if (bc == 'stress-free' and ell == 1):
        # Add null space solution
        λ_analytic = np.append(0, λ_analytic)
    assert np.allclose(solver.eigenvalues[:Nr//4], λ_analytic[:Nr//4])