DedalusProject / dedalus

A flexible framework for solving PDEs with modern spectral methods.
http://dedalus-project.org/
GNU General Public License v3.0
489 stars 115 forks source link

Incorrect conversion for cartesian NCC #281

Closed lecoanet closed 7 months ago

lecoanet commented 7 months ago

I'm trying to solve the following eigenvalue problem using Cartesian coordinates where u depends on y (ComplexFourier) and z (ChebyshevT) only:

import numpy as np
import dedalus.public as d3

ny = 2
nz = 384
dtype = np.complex128

coords = d3.CartesianCoordinates('x','y','z')
dist = d3.Distributor(coords, dtype=dtype)
ybasis = d3.ComplexFourier(coords['y'],size=ny,bounds=(0,1))
zbasis = d3.ChebyshevT(coords['z'],size=nz,bounds=(0.5,1))
basis = (ybasis,zbasis)

u = dist.VectorField(coords, bases=basis, name='u')
om = dist.Field(name='om')
v0 = dist.VectorField(coords, bases=zbasis, name='v0')
v0['g'][0] = 1

problem = d3.EVP([u], eigenvalue=om, namespace=locals())
problem.add_equation("om*u + cross(curl(u), v0) = 0")

solver = problem.build_solver()
subproblem = solver.subproblems_by_group[(0, 0, None)]
solver.solve_sparse(subproblem, 1, 0.)

print('solved eigenvalues',solver.eigenvalues)

When trying to build the matrices, I run into the following error:

  File "/Users/lecoanet/projects/dedalus/d3_new3/dedalus/core/basis.py", line 624, in _last_axis_component_ncc_matrix
    convert = jacobi.conversion_matrix(Nmat, arg_basis.a, arg_basis.b, out_basis.a, out_basis.b)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/lecoanet/projects/dedalus/d3_new3/dedalus/tools/jacobi.py", line 235, in conversion_matrix
    raise ValueError("a0 must be less than or equal to a1")
ValueError: a0 must be less than or equal to a1

This seems to be related to both the cross product and the curl operations.

Kyle Augustson found this issue.

lecoanet commented 7 months ago

@kburns pointed out that the error goes away if we change the basis of the NCC to zbasis.derivative_basis(1). Of course at that stage the matrix is singular, but that is the correct behavior. Nevertheless, the code should apply the correct conversion to the NCC without the user having to change the NCC basis.

kburns commented 7 months ago

This should be fixed with https://github.com/DedalusProject/dedalus/commit/b00863365cc230244dbe6e5e171445fb5b376367.

lecoanet commented 7 months ago

Kyle confirms this fixes the issue.