Possibly an extremely basic question about implementing the 2D vorticity evolution for the incompressive Navier-Stokes

[Adasfs Dsdade]

Hey all, I'm a non-domain-expert who wants to play around with some of this stuff, right now I'm seeing if I can implement the 2D vorticity evolution of the Navier-Stokes.

I've been having what is probably a very basic issue, but I can't seem to modify the example code for 2D shear flow to evolve through the vorticity equation rather than the Navier-Stokes. I keep getting linear operators on LHS type errors, but no matter how I shuffle the terms around I can't get this error to go away. Any hints?

Modified code snippet in question. Error is thrown for the bolded line.

# Fields
p = dist.Field(name='p', bases=(xbasis,zbasis))
s = dist.Field(name='s', bases=(xbasis,zbasis))
u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
w = dist.Field(name='w', bases=(xbasis,zbasis))
tau_p = dist.Field(name='tau_p')

# Substitutions
nu = 1 / Reynolds
D = nu / Schmidt
x, z = dist.local_grids(xbasis, zbasis)
ex, ez = coords.unit_vector_fields(dist)

# Problem
problem = d3.IVP([u, s, p, tau_p], namespace=locals())
problem.add_equation("div(skew(u)) = - w")
problem.add_equation("dt(w) = nu*lap(w) - u@grad(w)")
problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0") # Pressure gauge

[Jeffrey S. Oishi]

Hi,

Your problem variables are u, s, p, and tau_p. w is not a problem variable, thus dt(w) is not linear in the problem variables.

It would appear what you want is to define w as a substitution, rather than a problem variable: before the #problem section, do

w = d3.div(d3.skew(u))

Note also that your system is not square: you have five equations but only four variables. This will also fix that problem.

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/2a67ff74-6957-4363-b388-230e52c983ebn%40googlegroups.com.

[Adasfs Dsdade]

Hey Jeff, thanks for taking a look, I really appreciate it. So I updated the code per your suggestion (and tried a bunch of other things but I digress) to get

# Fields
p = dist.Field(name='p', bases=(xbasis,zbasis))
s = dist.Field(name='s', bases=(xbasis,zbasis))
u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))

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

# Substitutions
nu = 1 / Reynolds
D = nu / Schmidt
x, z = dist.local_grids(xbasis, zbasis)
ex, ez = coords.unit_vector_fields(dist)

w = -d3.div(d3.skew(u))

# Problem
problem = d3.IVP([u, s, p, tau_p], namespace=locals())

problem.add_equation("dt(w) - nu*lap(w) = - u@grad(w)")

problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0") # Pressure gauge

The output complained about a non-square system at the solver instantiation. OK, so as far as I can tell the dimensions on the LHS and RHS all match. I think the issue is the problem variables have one more dimension (u is 2D, whereas the matrix operators act on w which is 1D). When I make the fix

# Problem

problem = d3.IVP([w, s, p, tau_p], namespace=locals())
problem.add_equation("dt(w) - nu*lap(w) = - u@grad(w)")

problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")

problem.add_equation("tau_p = - div(u)")

problem.add_equation("integ(p) = 0") # Pressure gauge

I get a linear problem variable error. Since it's visibly linear (I think), I assume the problem is putting a substitution as a problem variable. With this in mind, I'm thinking maybe the move is to start with w instead of u, and recover u through the Biot-Savart law. As a last point, when I try

# Problem
problem = d3.IVP([u, s, p, tau_p], namespace=locals())

problem.add_equation("dt(w) - nu*lap(w) = - u@grad(w)")

problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")

problem.add_equation("w + div(skew(u)) = 0")

problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0") # Pressure gauge

I get a singularity error after a time step (expected since it just repeats the substitution), however this means it got past the square system problem! This is good!!! My takeaway is I need to add a line here that does nothing to the system without throwing singularities into the evolution. I'll experiment some more but I'd appreciate any more ideas you might have.

[Keaton Burns]

Hello,

This system is still not well posed for multiple reasons. Adding a constant to u will not change the LHS, because u only appears inside derivatives. Also no pressure modes appear in the linear system besides the mean mode, so they are unconstrainted. If you’re using a vorticity equation, you don’t need a pressure variable at all. You should drop p from the system, define w using a substitution, and add a gauge constraint on u.

Best,

-Keaton

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/660e11e1-5f56-498b-b530-192264ad0699n%40googlegroups.com.