Error with Curl in Shell Basis

[Oceanview24]

I have a simple script to test the curl operator in a shell basis:

# =============================================================================
# Minimal Failing Example for Dedalus v3.0.4 ShellBasis
#
# This script is designed to test the 'curl' operator. It is a mathematically
# well-posed problem (3 variables, 3 boundary conditions) but fails during
# solver build with a 'Non-square system' error, indicating a probable bug
# in the Dedalus library's handling of the curl operator in this geometry.
# =============================================================================
import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# --- Parameters ---
Nphi, Ntheta, Nr = 16, 8, 8 # Low resolution for a fast test
r_inner, r_outer = 1.0, 2.0
dtype = np.float64

# --- Setup ---
coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype)
shell = d3.ShellBasis(coords, shape=(Nphi, Ntheta, Nr), radii=(r_inner, r_outer), dealias=3/2, dtype=dtype)

# --- Fields ---
B = dist.VectorField(coords, name='B', bases=shell)

# --- Problem ---
problem = d3.IVP([B], namespace=locals())
# Add the simplest possible equation that uses the curl operator.
problem.add_equation("dt(B) + curl(B) = 0")
# Add the correct number of boundary conditions (3 for a 1st-order system of 3 variables).
problem.add_equation("B(r=r_outer) = 0")


# --- Build Solver ---
try:
    logger.info("Attempting to build the solver with a 'curl' operator...")
    solver = problem.build_solver(d3.RK222)
    logger.info("Solver built successfully.")
except Exception as e:
    logger.error("Solver build FAILED.")
    logger.error(f"Error Type: {type(e).__name__}")
    logger.error(f"Error Message: {e}")

However, I am getting this error even though the problem is well-posed:
2025-09-16 21:58:56,007 __main__ 0/1 INFO :: Attempting to build the solver with a 'curl' operator...
2025-09-16 21:58:56,036 __main__ 0/1 ERROR :: Solver build FAILED.
2025-09-16 21:58:56,036 __main__ 0/1 ERROR :: Error Type: ValueError
2025-09-16 21:58:56,036 __main__ 0/1 ERROR :: Error Message: Non-square system: group=(0, 0, None), I=9, J=8

This seems to be an issue with how Dedalus handles the curl operator in a shell basis.

[Jeffrey S. Oishi]

Hi,

This is not a bug. Your system is not well posed, as it does not specify the constant term of B. This has nothing to do with spherical coordinates; a simple test in Cartesian will fail in the same way. You need to add a tau variable. For more details, please see our documentation:

https://dedalus-project.readthedocs.io/en/latest/pages/tau_method.html

A corrected version of your script is attached. 

Jeff

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/04988252-e922-453c-97de-faf9acd1175an%40googlegroups.com.

[Jeffrey S. Oishi]

Oops, I should be more careful; your system is well-posed (the boundary condition specifies the constant), but the system is not square because you did not specify a variable to hold the boundary condition. 

[Oceanview24]

Hi,

I would like to add a condition that divergence of B = 0 everywhere, and add another boundary condition for the inner surface. I used the example given here:
https://dedalus-project.readthedocs.io/en/latest/pages/examples/ivp_shell_convection.html

I have attached my code, however it is still giving a non-square system I=18, J=11

How can I fix this?

Thanks!

[Jeffrey S. Oishi]

Hi,

You can't do this and have a square system unless you add an unphysical term to the induction equation. Incompressible hydro has a pressure term that enforces a divergence free condition on the velocity field. The induction equation has no equivalent term. The statement that the divergence of B is zero is an observational statement about the universe. While the induction equation will not introduce divergence of B if it is not present in the initial conditions, numerical approximations often will. There are many, many ways to deal with this (one of the simplest is to evolve the vector potential A instead; then you can enforce a gauge like div(A) = 0 with a pressure-like term that does not change the physics since the curl of a gradient is zero); you'll have to decide which is best for your particular purpose. 

Jeff

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/82239b12-a419-45ae-9127-cd05bed51bacn%40googlegroups.com.