Issue with EVP ("Objects are not all equal.")

[Christopher Li]

Hi,

I am trying to solve an EVP for linearized shallow water equation, but I kept getting error when adding the first equation to the problem. 

I noticed that if I do "dt = lambda A: omega*A", instead of "dt = lambda A: -1j*omega*A", the code runs. Why does this "i" stop the code from running?

Can someone help? Thanks!

Here are the code and the error message:

import numpy as np
import dedalus.public as d3
import logging
import matplotlib.pyplot as plt
import h5py
from dedalus.extras.plot_tools import plot_bot_2d
logger = logging.getLogger(__name__)

# Parameters
Lx, Ly = (10, 20*np.pi)
Nx, Ny = (128, 128)
dealias = 3/2
dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

# Fields
x, y = dist.local_grids(xbasis, ybasis)
ex, ey = coords.unit_vector_fields(dist)
u_vec = dist.VectorField(coords, name='u_vec', bases=(xbasis,ybasis))
u = dist.Field(name='u', bases=(xbasis,ybasis))
v = dist.Field(name='v', bases=(xbasis,ybasis))
eta = dist.Field(name='eta', bases=(xbasis,ybasis))
f = dist.Field(name='f', bases=(xbasis,ybasis))
omega = dist.Field(name='omega')
f['g'] = np.ones(Ny)*np.sin(2*np.pi*y/Ly)
# print('shape of f is',np.shape(f['g']))

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
# dx = lambda A: 1j*kx*A
# dy = lambda A: 1j*ky*A
dt = lambda A: -1j*omega*A

# Problem
problem = d3.EVP([u, v, eta], eigenvalue= omega, namespace=locals())
problem.add_equation("dt(u) + dx(eta) - f*v = 0")
problem.add_equation("dt(v) + dy(eta) + f*u = 0")
problem.add_equation("dt(eta) + dx(u) + dy(v) = 0")

-----------------------------------------------------------------------------------------------------

Error message:

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[29], line 43 41 # Problem 42 problem = d3.EVP([u, v, eta], eigenvalue= omega, namespace=locals()) ---> 43 problem.add_equation("dt(u) + dx(eta) - f*v = 0") 44 problem.add_equation("dt(v) + dy(eta) + f*u = 0") 45 problem.add_equation("dt(eta) + dx(u) + dy(v) = 0") File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/problems.py:75, in ProblemBase.add_equation(self, equation, condition) 73 namespace = dict(self.namespace) 74 LHS_str, RHS_str = parsing.split_equation(equation) ---> 75 LHS = eval(LHS_str, namespace) 76 RHS = eval(RHS_str, namespace) 77 else: 78 # Split operator tuples File <string>:1 File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/field.py:89, in Operand.__add__(self, other) 86 def __add__(self, other): 87 # Call: self + other 88 from .arithmetic import Add ---> 89 return Add(self, other) File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/tools/dispatch.py:17, in MultiClass.__call__(cls, *args, **kw) 14 """Dispatch instantiation based on the supplied arguments.""" 15 try: 16 # Preprocess arguments and keywords ---> 17 args, kw = cls._preprocess_args(*args, **kw) 18 # Create instance if no subclasses 19 subclasses = [sc for sc in cls.__subclasses__() if not getattr(sc, "stop_dispatch", False)] File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/arithmetic.py:69, in Add._preprocess_args(cls, *args, **kw) 67 dist = unify_attributes(args, 'dist', require=False) 68 tensorsig = unify_attributes(args, 'tensorsig', require=False) ---> 69 dtype = unify_attributes(args, 'dtype', require=False) 70 args = [Operand.cast(arg, dist, tensorsig, dtype) for arg in args] 71 return args, kw File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/tools/general.py:85, in unify_attributes(objects, attr, require) 83 if require: 84 raise ---> 85 return unify(attrs) File ~/miniconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/tools/general.py:72, in unify(objects) 70 else: 71 if object != OBJECT: ---> 72 raise ValueError("Objects are not all equal.") 73 return OBJECT ValueError: Objects are not all equal.

[Ben Brown]

Chistopher,

     You need to set dtype = np.complex128 rather than dtype = np.float64.  You can see examples of this in e.g., evp_1d_mathieu/mathieu_evp.py in our examples directory.

The issue is from introducing a complex omega on a real grid.

--Ben

--
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/ddab5b5d-a931-4197-82a0-eebcc40d16c1n%40googlegroups.com.