complicated boundary conditions

[Alexander Morozov]

Dear all

We have recently started using Dedalus to simulate viscoelastic flows. Thanks a lot for a great piece of software!

I have a conceptual question. We are simulating some 2D flow on a Fourier-Chebyshev domain. Symbolically, the equation of motion looks like this

d_t f = Laplacian(f) + nonlinear terms(f)                  (*)

where f's are a bunch of 2D fields. Obviousy, it requires two boundary conditions for each field f. However, these are unknown. What is typically done in the literature, is to evaluate (*) in the bulk, while solving the same equation without the Laplacian at the boundaries, i.e.

d_t fb = nonlinear terms(f)        (at each boundary)               (**) 

Here fb only lives in 1D (along the Fourier direction). Then one can use the solutions too (**) as the boundary conditions for (*).

Do you have any suggestions as to how this can be implemented in Dedalus?

I was thinking of creating two 1D fields, say, f_up and f_down, updating them according to (**), and feeding them to the original problem as bc. This is not particularly elegant. Can this be done? Or would you have a better idea?

Thank you and best wishes,

Alexander

[Alexander Morozov]

I guess my question boils down to this: can one simultaneously have 2D and 1D fields defined on the same (2D) domain?

Thanks,

Alexander 

[Keaton Burns]

Hi Alexander,

Yes it’s possible to add fields and set the problem metadata to indicate that they are constant in the z-direction, and then solve 1D equations for them.  However at first glance that seems like it might be a poorly conditioned formulation. If the boundary condition just becomes “f = fb” at the boundary, then differentiating this equation in time and subtracting the evolution equations would indicate that this is the same boundary condition as imposing “Laplacian(f) = 0” at the boundaries, right?  And I don’t think that works well as a BC for the Laplacian itself.  Or are there other terms from the RHS that survive?  Just curious and don’t have much experience with viscoelastic simulations, myself!

Best,

-Keaton

--
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/a0aa8856-8849-4a0e-b650-6400d889cb2en%40googlegroups.com.

[Evan Henry Anders]

Hi Alexander,

Yes, using interpolation you can get horizontally 1D fields (f_up and f_down) from f at the boundaries.  Say you have an equation like:

problem.add_equation("dt(f) - D*(dx(dx(f)) + dz(fz))       = -(u*dx(f) + w*fz)")

I believe you can straightforwardly add boundary conditions like:

problem.add_bc("left(dt(b)) = left(-(u*dx(b) + w*bz))")

problem.add_bc("right(dt(b)) = right(-(u*dx(b) + w*bz))")

Where the "left" equations is f_down at the lower boundar and the "right" equations is f_up at the upper boundary. This does what I *think* you want to do. But, like Keaton, I have very little experience with viscoelastic simulations. I think the above will timestep, but I don't know if its results will make sense! Curious to hear a bit more about your problem, and hopefully this helps you get started towards a solution.

Best,

Evan

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBiJa6gbbCkdKitoYaFMC6%3DdBphFoL0NFW2F5VtZa7AZgQ%40mail.gmail.com.

[Evan Henry Anders]

Hi Alexander,

Apologies, in my "add_bc" equations, I meant to replace "b" with "f" but I didn't proofread my email very carefully. Those should read:

problem.add_bc("left(dt(f)) = left(-(u*dx(f) + w*fz))")

problem.add_bc("right(dt(f)) = right(-(u*dx(f) + w*fz))")

Where u and w are some horizontal and vertical velocity fields and fz = dz(f).

Best,

Evan

[Alexander Morozov]

Hi Keaton & Evan

Thank you for your responses. I currently do the following:

# Create bases and domain

x_basis = de.Fourier('x', nx, interval=(0, Lx), dealias=3/2)

y_basis = de.Chebyshev('y',ny, interval=(-1, 1), dealias=3/2)

domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)

problem = de.IVP(domain, variables=['p','u','v','uy','vy','cxx', 'cxy', 'cyy', 'Dcxx', 'Dcxy', 'Dcyy', 'gxx', 'gxy', 'gyy'])

problem.meta[:]['y']['dirichlet'] = True

problem.parameters['Remin1'] = 1.0/Reynolds

problem.parameters['Wimin1'] = 1.0/Wi

problem.parameters['beta'] = beta

problem.parameters['Lx'] = Lx

problem.parameters['Ly'] = Ly

problem.parameters['D'] = diffusivity

problem.parameters['FENEL2'] = FENEL*FENEL

problem.substitutions['TrC'] = "1.0 - (cxx+cyy)/FENEL2"

problem.add_equation("dt(u) + dx(p) - beta*Remin1*(dx(dx(u)) + dy(uy)) \

        = Remin1*Wimin1*(1.0-beta)*( dx(cxx/TrC) + dy(cxy/TrC) ) - u*dx(u) - v*uy + 2*Remin1")

problem.add_equation("dt(v) + dy(p) - beta*Remin1*(dx(dx(v)) + dy(vy)) \

        = Remin1*Wimin1*(1.0-beta)*( dx(cxy/TrC) + dy(cyy/TrC) ) - u*dx(v) - v*vy")

problem.add_equation("dx(u) + vy = 0")

problem.add_equation("dt(cxx) - D * ( dx(dx(cxx))+dy(Dcxx) ) = - u*dx(cxx) - v*Dcxx + 2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1")

problem.add_equation("dt(cxy) - D * ( dx(dx(cxy))+dy(Dcxy) ) = - u*dx(cxy) - v*Dcxy + cxx*dx(v) + cyy*uy - Wimin1*cxy/TrC")

problem.add_equation("dt(cyy) - D * ( dx(dx(cyy))+dy(Dcyy) ) = - u*dx(cyy) - v*Dcyy + 2.0*cxy*dx(v) + 2.0*cyy*vy - Wimin1*cyy/TrC + Wimin1")

problem.add_equation("dt(gxx) = - u*dx(cxx) - v*Dcxx + 2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1")

problem.add_equation("dt(gxy) = - u*dx(cxy) - v*Dcxy + cxx*dx(v) + cyy*uy - Wimin1*cxy/TrC")

problem.add_equation("dt(gyy) = - u*dx(cyy) - v*Dcyy + 2.0*cxy*dx(v) + 2.0*cyy*vy - Wimin1*cyy/TrC + Wimin1")

problem.add_equation("uy - dy(u) = 0")

problem.add_equation("vy - dy(v) = 0")

problem.add_equation("Dcxx - dy(cxx) = 0")

problem.add_equation("Dcxy - dy(cxy) = 0")

problem.add_equation("Dcyy - dy(cyy) = 0")

problem.add_bc("left(u) = 0")

problem.add_bc("right(u) = 0")

problem.add_bc("left(v) = 0")

problem.add_bc("right(v) = 0", condition="(nx != 0)")

problem.add_bc("left(p) = 0", condition="(nx == 0)")

problem.add_bc("left(cxx) - left(gxx) = 0")

problem.add_bc("right(cxx) - right(gxx) = 0")

problem.add_bc("left(cxy) - left(gxy) = 0")

problem.add_bc("right(cxy) - right(gxy) = 0")

problem.add_bc("left(cyy) - left(gyy) = 0")

problem.add_bc("right(cyy) - right(gyy) = 0")

This is not particularly optimal as it increases the problem size by the additional dummy fields gxx, gxy, and gyy that I only need to know at the boundary. 

[Alexander Morozov]

Hi Keaton

What one does is a bit more subtle than Laplacian(f) = 0. It's easier to understand in terms of operator splitting. In our own code,  at the wall, we first timestep all the terms without the Laplacian to some intermediate time t* (a fraction of the full timestep). Then you timestep only the Laplacian term to the full timestep. This needs two bc and one uses the boundary values from the intermediate solution at t* there. This allows one to simultaneously include the stabilising Laplacian in the code and removes the need to invent new b.c. that no-one really knows what they should be. 

Alexander

[Alexander Morozov]

Hi Evan

I am very new to this, can I please check that I understand this properly. I'll use my particular example:

problem.add_equation("dt(cxx) - D * ( dx(dx(cxx))+dy(Dcxx) ) = - u*dx(cxx) - v*Dcxx + 2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1")

Can I just say 

problem.add_bc("left(dt(cxx))  = left(2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1)")

and similar for other fields/boundaries?

Best wishes,

Alexander

[Evan Henry Anders]

Hi Alexander,

Yep! Exactly as you've written it is how I would try setting up this problem myself.

One (completely optional) thing that you can do to ensure that your equations are consistent is that you can use substitutions, e.g.,

problem.substitutions['cxx_forcing'] = '2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1'
problem.add_equation("dt(cxx) - D * ( dx(dx(cxx))+dy(Dcxx) ) = - u*dx(cxx) - v*Dcxx + cxx_forcing")
problem.add_bc("left(dt(cxx))  = left(cxx_forcing)")

When specifying long strings that are repeated in two places, substitutions like this can be a great way to make sure that there are no typo'd differences between the string in the two equations.  

Best,

Evan

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/87e84e13-07bd-4613-b946-609a525c1065n%40googlegroups.com.

[Alexander Morozov]

Hi Evan

Thanks again for your help. Here is a short update. I have spent some time studying the following two options:

problem.substitutions['cxx_forcing'] = '2.0*cxx*dx(u) + 2.0*cxy*uy - Wimin1*cxx/TrC + Wimin1'

Option 1: my original solution

problem.add_equation("dt(cxx) - D * ( dx(dx(cxx))+dy(Dcxx) ) = - u*dx(cxx) - v*Dcxx + cxx_forcing")

problem.add_equation("dt(gxx) =  cxx_forcing")

problem.add_bc("left(cxx) - left(gxx) = 0")

etc...

Option 2: following your suggestion

problem.add_equation("dt(cxx) - D * ( dx(dx(cxx))+dy(Dcxx) ) = - u*dx(cxx) - v*Dcxx + cxx_forcing")
problem.add_bc("left(dt(cxx))  = left(cxx_forcing)")

etc...

Option 1 is rather ugly as it introduces an extra field g for every physical field c. As g's are only needed at the boundaries, this feels excessive  and rather expensive, while Option 2 looks elegant. However, when I compared the code performance of these options, it appears that Option 1 is about 4 times faster than Option 2. These have been obtained for 256x1024 and 512x2048 Fourier/Chebyshev resolutions on 16/32/64 cores on a single node (2X AMD 7452).

I guess the interpolation towards the boundary does not result in a nice linear operator(?) but a factor of ~four is surprising. I would have certainly guessed the other way round. 

Thanks again and best wishes,

Alexander

[Keaton Burns]

Hi Alexander,

Hmm that is definitely unexpected.  If you’d like to share some minimal working scripts showing this difference, I can do a little profiling to see what’s causing the slowdown.

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/36131efa-3396-4f85-a3c6-c5c8df64e40cn%40googlegroups.com.

[Alexander Morozov]

Hi Keaton

Please see attached. The file names correspond to the options given in my previous post. Run times on 32 cores: option 1 - 15 sec, option 2 - 119 s.

BTW does it matter whether one does left(A+B) = left(C*D) or left(A) + left(B) = left(C)*left(D)?

Thank you,

Alexander

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBhYOytfN6t5hBeZ%2B9kN3gb2YJThSRu_y%2B0c0%3DWr6odmzw%40mail.gmail.com.