Issue with Non-constant Coefficients in an EVP

[Nathan Stolnicki]

Hello!

I have a quick question regarding non-constant coefficients in an EVP. I am attempting to solve the following eigenvalue problem:

coords = d3.CartesianCoordinates('x','y','z')

dist = d3.Distributor(coords, dtype=np.complex128)

xbasis = d3.ComplexFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=3/2)

ybasis = d3.ComplexFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=3/2)

x, y = dist.local_grids(xbasis, ybasis)

n = dist.Field(name='n', bases=(xbasis,ybasis))

v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis))

B = dist.VectorField(coords, name='B', bases=(xbasis,ybasis))

B_p = dist.VectorField(coords, name='B_p', bases=(xbasis,ybasis))

J = dist.VectorField(coords, name='J', bases=(xbasis,ybasis))

omega = dist.Field(name='omega')

dx = lambda A: d3.Differentiate(A, coords['x'])

dy = lambda A: 1j*ky*A

dt = lambda A: -1j*omega*A

B_p = B - d_e*d_e*d3.Laplacian(B)

J = d3.Curl(B)

B0 = dist.VectorField(coords,name='B0',bases=(xbasis,ybasis))

B0['g'][0] = np.sin(2*np.pi/Lx * y)

n0 = dist.Field(name='n0',bases=(xbasis,ybasis))

n0['g'] = -1/2 * np.sin(2*np.pi/Lx * y)*np.sin(2*np.pi/Lx * y) + 1

problem = d3.EigenvalueProblem(variables=[B, n, v], eigenvalue=omega, namespace=locals())

problem.add_equation("dt(n) + n0*div(v) = 0")

problem.add_equation("n0*dt(v) + grad(T*n + B0@B) + B0@grad(B) = 0")

problem.add_equation("dt(B_p) - B0@grad(v) + B0*div(v) - (1/n0)*B0@grad(J) = 0")

However, I obtain the following error:

 “dedalus.tools.exceptions.SymbolicParsingError: Must build NCC matrices with same variables”.

This error persists even if the nccs, B0 and n0, are not set and remain at zero. Thus, I presume the issue lies in how the nccs (especially B0) are initialized. I have never used a ncc vector field in cartesian coordinates before in Dedalus so any help on my implementation would be greatly appreciated! 

Best,

Nathan

[Daniel Lecoanet]

Hi Nathan,

I cannot reproduce your error because you did not include the entire script (e.g., what are Nx, Ny, T?). However, I note that B0 and n0 are both functions of only y, so you should define them as:

B0 = dist.VectorField(coords,name='B0',bases=ybasis)

n0 = dist.Field(name='n0',bases=ybasis)

If you attach the entire code I could try to reproduce the error.

Daniel

--
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 on the web visit https://groups.google.com/d/msgid/dedalus-users/1325bee8-f89d-4a20-906b-e1f18e1a34e1n%40googlegroups.com.

[Nathan Stolnicki]

Hi Daniel,

Thank you for the quick response! Attached is the full version of my code. I tried your suggestion about changing B0 and n0 to only have the ybasis, however, I still recieve the same error message as above. Thank you for any help you can provide!

Best,

Nathan

[Daniel Lecoanet]

Hi Nathan,

Thanks, I’m able to reproduce the issue, and I believe it might be a bug in Dedalus. We will investigate further and get back to you.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/847c5cc3-c144-4bfa-a407-94a89793ef85n%40googlegroups.com.
<EVP.py>

[Daniel Lecoanet]

Hi Nathan,

As I noted, there is a bug associated with taking the gradient of the curl. A temporary fix is to switch your equation for B_p to this:

problem.add_equation("dt(B_p) - B0@grad(v) + B0*div(v) + (1/n0)*B0@grad(-J) = 0”)

If you have a negative sign inside the gradient, the code will do the proper thing. So I would use that workaround for now.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CDAB8B41-6989-40EF-BDA3-B0F433930BDA%40northwestern.edu.