[d3] Terms on RHS not always properly type-cast to same dtype as rest of domain

[bpbrown]

I've found an interesting edge case where something in d3 didn't behave as I expected. I'm not sure if this is a valid edge case that we need to address, or just a spot where my expectations were wrong.

As background, I've been running an NLBVP problem on a complex domain that has a numpy ufunc on the RHS. The problem looks like this:

        dtype=np.complex128
        <...>
        b = basis.BallBasis(c, (1, 1, Nr), radius=1, dtype=dtype, dealias=dealias)
        <...>
        problem = problems.NLBVP([f, R, τ], ncc_cutoff=ncc_cutoff)
        problem.add_equation((lap(f) + LiftTau(τ), - R**2 * Pow(f,n)))
        problem.add_equation((f(r=0), 1))
        problem.add_equation((f(r=1), np.exp(-3/n)))

This is a truncated version of Lane-Emden that will capture the inner 3 density scale heights of the solution, and can be achieved by modifying the last boundary condition in the test_nlbvp.py unit test test_lane_emden_floating_R() .

This problem works fine when we use problem.add_equation((f(r=1), 0) but with problem.add_equation((f(r=1), np.exp(-3/n)) we get an interesting error:

  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/core/problems.py", line 160, in add_equation
    expr = LHS - RHS
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/core/field.py", line 67, in __sub__
    return (self + (-other))
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/core/field.py", line 58, in __add__
    return Add(self, other)
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/tools/dispatch.py", line 26, in __call__
    args, kw = cls._preprocess_args(*args, **kw)
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/core/arithmetic.py", line 52, in _preprocess_args
    dtype = unify_attributes(args, 'dtype', require=False)
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/tools/general.py", line 85, in unify_attributes
    return unify(attrs)
  File "/Users/bpbrown/conda_install/src/dedalus-d3/dedalus/tools/general.py", line 72, in unify
    raise ValueError("Objects are not all equal.")
ValueError: Objects are not all equal.

This ValueError can be fixed by doing either of two things. You can directly type-cast the ufunc:

        problem.add_equation((f(r=1), np.exp(-3/n, dtype=dtype)))

or you can define a correctly-typed Field object to store the value ahead of time:

    f_top = field.Field(dist=d, dtype=dtype, name='f_top')
    f_top['g'] = np.exp(-3/n) # this avoids problems with complex vs float
    problem.add_equation((f(r=1), f_top))

It seems like it would be nice if the RHS got automatically type-cast to the right dtype as the rest of the domain: the RHS seems like the most likely place where this edge case might come up, from constant values or boundary forcing.

But maybe doing this causes some broader set of problems, since fields are typed independently of the problem.

What are our thoughts?

[kburns]

Good catch. I think you're right that we want to just cast up here, but we have to be careful since we probably want it to still fail on the opposite case (the LHS is real operator expression and the RHS is given as a complex number).