Semi-infinite domain; time-dependent boundary conditions

[Richard Katz]

Hi all,

I'm a new Dedalus user seeking advice to get started on an IVP application code.  These are the points on which I'd be grateful for your guidance.

1 - Whether it is possible (and if so, how) to set up a semi-infinite, two dimensional Cartesian domain.  Periodic in x, semi-infinite in y.  

2 - How to impose boundary conditions at y=0 that depend explicitly time and position (x).

Thanks in advance for any help,

Richard

[Jeffrey S. Oishi]

Hi Richard,

1 - Whether it is possible (and if so, how) to set up a semi-infinite, two dimensional Cartesian domain.  Periodic in x, semi-infinite in y.  

Not currently. We have not implemented Laguerre polynomials in the current version of Dedalus (though they exist in version 2), though this is something we may do in the medium term.

 

2 - How to impose boundary conditions at y=0 that depend explicitly time and position (x).

That is straightforward. You can see an example here:

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

Despite it being a boundary value problem, the basic idea about boundary conditions remains the same in an initial value problem. 

Jeff

[Richard Katz]

Hi Jeff,

Thank you for your reply and pointing me to the helpful example.  I have made some progress but am getting an error.  Here’s a code snippet:

# Bases

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

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

xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias)

zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias)

# Fields

t = dist.Field()

p = dist.Field(name='p', bases=(xbasis,zbasis))

u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))

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

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

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

# Substitutions

x, z = dist.local_grids(xbasis, zbasis)

ex, ez = coords.unit_vector_fields(dist)

lift_basis = zbasis.derivative_basis(1)

lift = lambda A: d3.Lift(A, lift_basis, -1)

grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction

# Boundary condition

k = 2*np.pi/wavelength

omega = 2*np.pi/period

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

#U['g'][0]=np.cos(k*x)

U['g'][0]=np.cos(k*x-omega*t) # THIS IS THE KEY LINE

U['g'][1]=0

# Problem

# First-order form: "div(f)" becomes "trace(grad_f)"

# First-order form: "lap(f)" becomes "div(grad_f)"

problem = d3.IVP([p, u, tau_p, tau_u1, tau_u2], time=t, namespace=locals())

problem.add_equation("trace(grad_u) + tau_p = 0")

problem.add_equation("dt(u) - Nu*div(grad_u) + grad(p) + lift(tau_u2) = - u@grad(u)")

problem.add_equation("u(z=0) = U")

problem.add_equation("u(z=Lz) = 0")

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

And here’s the error:

Traceback (most recent call last):

  File "/Users/katz/Desktop/dedalus/travelling_stokes/travelling_stokes.py", line 52, in <module>

    U['g'][0]=np.cos(k*x-omega*t)

                     ~~~^~~~~~~~

  File "/Users/katz/Applications/miniforge3/envs/dedalus3/lib/python3.13/site-packages/dedalus/core/field.py", line 60, in __array_ufunc__

    return arithmetic.Add(inputs[0], (-1)*inputs[1])

           ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^

  File "/Users/katz/Applications/miniforge3/envs/dedalus3/lib/python3.13/site-packages/dedalus/tools/dispatch.py", line 17, in __call__

    args, kw = cls._preprocess_args(*args, **kw)

               ~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^

  File "/Users/katz/Applications/miniforge3/envs/dedalus3/lib/python3.13/site-packages/dedalus/core/arithmetic.py", line 58, in _preprocess_args

    args = [arg for arg in args if arg != 0]

                                   ^^^^^^^^

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

If I just use the boundary condition

U['g'][0]=np.cos(k*x)

Then the code runs (I haven’t examined the output…)

So it seems that there’s a problem combining two Field variables in the same equation… ?

Sorry if I’m missing something obvious.

Richard

--
You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/M2ujC8OrhMY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/CAAOZTN9wa1AKQYwuMjQQVBTB14iZoaT1vsqjudHyetdDBgnMVA%40mail.gmail.com.

[Daniel Lecoanet]

Hi Richard,

The issue is that t is a field, whereas U[‘g’][0] is a numpy array. You cannot set the value of a numpy array to a field. What would be better would be for you to define U to be an operator, which evaluates to cos(k x-om t) for any value of the field t. One way to do this would be to define

U = np.cos(k*x_f - omega*t)*ex

where x_f is a field equal to x. Then when t is updated, U will be evaluated to different values. Please look at the documentation if you have questions about operators vs fields.

That said, it is actually best not to define a field x_f where x_f[‘g’]=x because x_f is not periodic, so if any transform occur it will incur substantial transform error. So an even better approach would be to define fields coskx and sinkx, which are periodic, and then

U = (coskx*np.cos(omega*t) + sinkx*np.sin(omega*t))*ex

Hope that helps,

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 visit https://groups.google.com/d/msgid/dedalus-users/EF0CA2BD-15FC-409E-B1C5-EF12AC8E6271%40gmail.com.

[Richard Katz]

Hi Daniel,

Just to confirm,  I think what you’ve suggested is 

coskx = dist.Field(bases=xbasis)

sinkx = dist.Field(bases=xbasis)

coskx['g'] = np.cos(k*x)

sinkx['g'] = np.sin(k*x)

U = (coskx*np.cos(omega*t) + sinkx*np.sin(omega*t))*ex

which seems to give me what I wanted.

Please correct me if this is wrong.

Thank you for your help!

Richard

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/3C4C447F-C404-4801-A8E0-5A910F98DDC5%40northwestern.edu.

[Daniel Lecoanet]

Yes, that looks good. It’ll be even more efficient if you do

U =  (d3.Grid(coskx)*np.cos(omega*t) + d3.Grid(sinkx)*np.sin(omega*t))*ex

which prevents superfluous transforms of coskx and sinkx between coefficient and grid space.

Daniel

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/92BF9717-FE4D-4804-8A98-A13D35B5905A%40gmail.com.