IVP: time dependent boundary conditions & proper syntax for performing an x integral at specific z

[Drummond Fielding]

Hello all,

Yesterday Daniel sat down with me and helped me run my first real simulation with dedalus! We designed a problem to simulate a 2D plume in which matter flows in a narrow stream into an stratified medium. We ran into one issue plus I have one question which is probably trivial. First the problem, in order to avoid a discontinuity during the first time steps when the plume first enters the domain we tried to impose a time dependent bc that essentially ramped up the plume using a tanh smoothing. To be exact this is what we did:

problem.add_bc("left(b) = B*0.5*(tanh((x+a/2)/sigma)-tanh((x-a/2)/sigma))*0.5*(1-tanh((t0-t)/delt))" )

problem.add_bc("left(w) = W*0.5*(tanh((x+a/2)/sigma)-tanh((x-a/2)/sigma))*0.5*(1-tanh((t0-t)/delt))" )

The part in bold is temporal ramp up of the plume in the vertical velocity and buoyancy. Luckily for us the code actually worked without the smoothing so I was able to run my first simulation. Nevertheless I would like to know if this is possible to do properly. When the temporal smoothing is in place I got the following error: http://pastebin.com/Y1uyK2fv

In case you want to see the whole code here is the script we came up with yesterday (based on the 3D Rayleigh-Benard example): http://pastebin.com/dkU6M1Wd . The relevant lines are numbers 46 and 48. 

OK, so that is the main issue we had and now a presumably less complicated question. I would like to have a bc at the top of the domain so that the mass flux out of the top is equal to what is inject in the bottom. To do so I need to perform a horizontal integration at z = 0. My question is what is the proper syntax to do this. Here is the integral I am trying to input:

    int_x [ b(z=0)*w(z=0) ]  / int_x [w(z=0)]

Thanks in advance for you help, and thanks for developing such an awesome code!

best,

Drummond Fielding

[Keaton Burns]

Hi Drummond,

Your code for the temporal smoothing is right, and the error you ran into was due to a bug in some of my recently-pushed code.  I’ve added a patch to the main repository, so it should work now.  Thanks for the catch!

The syntax for that expression looks like this:

integ(interp(b*w, z=0), ‘x’) / integ(interp(w, z=0), ‘x’)

If you’d like, you can use the new substitutions feature to define macros for simplifying this a bit, i.e.

problem.substitutions[“int_x(var)"] = “integ(var, ‘x’)"

The syntax of using “w(z=0)” to generate interpolation is on my to-do list, so that should hopefully be available soon, as well.

Also I'd definitely be interested in hearing about how things go with the outflow boundary conditions — not something I’ve played around with too much, myself.

Thanks!

-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 post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/659652a8-ea8d-4d8e-beaa-0fd871976713%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Drummond Fielding]

thanks a lot keaton i will try this now

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/RZ_6gX-WQSI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.

To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/99A5447F-F1AF-4E92-A2F0-9C706C657E12%40gmail.com.

[Aaron W]

Hello all,

   I am trying to do something very similar to what Drummond has done with BCs, but in addition to functions of time, I am trying to modulate by some discrete function, e.g.

                   problem.add_bc("right(u) = - a * omega * sqrt(n**2 - omega**2)/omega * right(shape * cos(omega*t - kx*x))")  

where shape is a field added as a parameter (but really only depends on x).

   But now that shape is required to be a field, must I define the parity to decompose in the x direction? I receive the following error, which I think is related to the missing parity definition, but it's not evident how to add parity to this field/parameter.

Traceback (most recent call last):

  File "boussinesqchebcos2D.py", line 143, in <module>

    problem.add_bc("right(u) = - a * omega * sqrt(n**2 - omega**2)/omega * right(shape * cos(omega*t - kx*x))")

  File "/Users/awink/dedalus/dedalus/core/problems.py", line 125, in add_bc

    self._build_object_forms(temp)

  File "/Users/awink/dedalus/dedalus/core/problems.py", line 142, in _build_object_forms

    temp['RHS'] = future.FutureField.parse(temp['raw_RHS'], self.namespace, self.domain)

  File "/Users/awink/dedalus/dedalus/core/future.py", line 223, in parse

    expression = eval(string, namespace)

  File "<string>", line 1, in <module>

  File "/Users/awink/dedalus/dedalus/core/operators.py", line 1335, in right

    return basis.Interpolate(arg0, 'right', out=out)

  File "/Users/awink/dedalus/dedalus/core/operators.py", line 1257, in __new__

    elif (arg0.meta[cls.basis.name]['constant']):

  File "/Users/awink/dedalus/dedalus/tools/cache.py", line 29, in __get__

    attribute = self.method(instance)

  File "/Users/awink/dedalus/dedalus/core/future.py", line 178, in meta

    meta[axis][key] = getattr(self, 'meta_%s' %key)(axis)

  File "/Users/awink/dedalus/dedalus/core/operators.py", line 649, in meta_parity

    return parity0 * parity1

TypeError: unsupported operand type(s) for *: 'NoneType' and 'int'

(I had previously accomplished this with the hacky helper function definitions defined after the problem definition.)

My new script for dedalus(3) is included below:

Thank you!

   Aaron

start_init_time = time.time()

x_basis = de.SinCos('x', Nx, interval=(-lx/2, lx/2), dealias=3/2)

z_basis = de.Chebyshev('z',Nz, interval=(-lz/2, lz/2), dealias=3/2)

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

x = domain.grid(0)

z = domain.grid(1)

xg, zg = np.meshgrid(x,np.transpose(z))

envelope = shapes.envelope(x,'Tophat',width,position).T

shapefield = np.zeros_like(xg)

for i in range(Nz):

    shapefield[i,:] = envelope

shape = domain.new_field()

shape['g'] = shapefield.T

# 2D Boussinesq Hydrodynamics

problem = de.IVP(domain, variables=['u','w','P','rho','u_z','w_z','rho_z'],time='t')

problem.meta['P','w','w_z','rho','rho_z']['x']['parity'] = 1

problem.meta['u','u_z']['x']['parity'] = -1

problem.parameters['Re'] = Re          # = \lambda*v_phase/\nu

problem.parameters['Fr'] = Fr           # = v_phase^2/(g*l\ambda)

problem.parameters['Sc'] = Sc           # = \nu/D

problem.parameters['eta'] = eta         # d\rho/dz

problem.parameters['omega'] = omega     # = 2*pi

problem.parameters['n'] = n             # = N/Omega

problem.parameters['shape'] = shape

problem.parameters['a'] = a

problem.parameters['kx'] = kx

problem.add_equation("dt(u)   - 1/Re*(dx(dx(u)) + dz(u_z)) + dx(P)               = - u*dx(u)   - w*dz(u)")

problem.add_equation("dt(w)   - 1/Re*(dx(dx(w)) + dz(w_z)) + dz(P)  + 1/Fr*(rho) = - u*dx(w)   - w*dz(w)") # Momentum

problem.add_equation("dt(rho) - 1/(Sc*Re)*(dx(dx(rho)) + dz(rho_z)) + w*eta      = - u*dx(rho) - w*dz(rho)") # Density

problem.add_equation("dx(u) + w_z = 0")  # Incompressible Continuity

problem.add_equation("dz(w)   - w_z   = 0")

problem.add_equation("dz(u)   - u_z   = 0")

problem.add_equation("dz(rho) - rho_z = 0") # Derivative Definitions

problem.add_bc("left(u) = 0") # Stress free: u_z = 0, dx(w) = 0

problem.add_bc("left(w) = 0", condition="(dx != 0)") # Need 6 BCs in z for each of 6 z differential equations

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

problem.add_bc("right(u) = - a * omega * sqrt(n**2 - omega**2)/omega * right(shape * cos(omega*t - kx*x))")

problem.add_bc("right(w) = a * omega * right(shape * cos(omega*t - kx*x))")

problem.add_bc("right(rho) = - a * eta * right(shape * sin(omega*t - kx*x))")

problem.add_bc("integ_z(P) = 0", condition="(dx == 0)") # Fixed Boundaries (Uses freed up mode, dx=i*k_x=0, to set gauge P.)

solver = problem.build_solver(de.timesteppers.SBDF3)

[Daniel Lecoanet]

Hi Aaron,

Yes, I think the missing thing is to define the parity of shape.  You can do this by setting problem.meta['shape']['x']['parity']=+-1.  I think this needs to be done after you define the parameter shape.  Also, I wonder if your BC's satisfy the parity constraint -- rho & w are supposed to be representable with only sin(x) modes, which is not going to be possible if they are equal to sin(kx x - omega t) on one boundary.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/4c08ac63-4adb-4b8d-af6c-651d56f87c6e%40googlegroups.com.

[Aaron W]

Hi Daniel,

   I appreciate your quick reply.

I had tried that, to which I get a key error since I don't think the parameters are being treated in the same way as the fields are as far as initialising the metadata structure. I thought a temporary workaround would work with problem.__setitem__()  even though that's horrible form, but I'll see where that takes me...

My BCs should satisfy parity only because of the parameter 'shape' (being locally confined and -> 0 towards the cos/sin boundaries). So I don't think this should be a problem---at least if I understood your point correctly.

Thank you!

   Aaron

[Daniel Lecoanet]

Ah, ok, I think I didn't say correctly how to set the metadata.  I think the following should work:

shape = domain.new_field()

shape['g'] = shapefield.T

shape.meta['x']['parity'] = -1

...

problem.parameters['shape'] = shape

Let me know whether or not that works.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/2677a2cd-e1b7-4116-9b9c-59ab43b1d6d0%40googlegroups.com.

[Keaton Burns]

Hi Aaron,

I think you might have to define two shape fields, one with each parity, to use in the boundary conditions for the fields with each respective parity.  The wave-like expressions should both be interpreted as having positive parity (since this is how ‘x' is interpreted), so I think the envelope for ‘u' would need to have negative parity and the envelopes for ‘w' and ‘rho' would need to have positive parity for the all of the BCs to match the fields.

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAJoYf%3DjC%2BzjW8ZA9_mjvPvZ_%2Bpd6Ujd3%3DUGsScUWWeF3oJmazg%40mail.gmail.com.

[Aaron W]

Daniel and Keaton,

   Thank you very much, Daniel -- using shape.meta['x']['parity']  worked perfectly!

And that is a really good point Keaton (and probable source of my parity mismatch errors)! So in summary, the base parity of the field should match the parity of it's respective BC(s)

Thanks!

   Aaron