Adding forcing to 2-D fluid equations

[Jonathan Squire]

Hi,

I'm a grad student at PPPL in Princeton, and have been starting to use dedalus for a little 2-D fluids side project. Geoff and I had a really interesting chat a few days ago in Sydney and it seems like you have a great project going!

My problem is 2-D the Charney-Hasegawa-Mima equation (or quasi-geostrophic equation) with hard wall conditions in one dimension, perfect for dedalus. However, I couldn't figure out how to add forcing into the system. The exact form of the forcing is not too important, but ideally I would be able to use both white noise (with some chosen spatial correlation), or a prescribed deterministic forcing. Geoff mentioned he may have a student looking into a nice implementation of such a feature soon, but is there a simple way to add something like this in the meantime? 

I've attached a basic script (mostly copied from your convection example). The script works fine, but doesn't have the forcing in it (there are some vain attempts to use evaluators to achieve this that are commented out). Any help would be much appreciated!

Thanks a lot,

Jono

PS. I had one other problem, related to using periodic boundary conditions in both directions. Basically, if I set the x_basis = de.Chebyshev(... to  x_basis = Fourier(... and remove the boundary conditions it gives me errors about singular matrices. Is there some way to use 2-D periodic boundary conditions (mainly for validation that my script is working correctly)?

[Ben Brown]

Jonathan,

    Welcome!

To answer your ps question: yes, you can do doubly periodic, but if you have a gauge condition you'll need to add an internal BC.

There's an example here:

https://bitbucket.org/bpbrown/kelvin_helmholtz

See especially KH_periodic.py and

https://bitbucket.org/bpbrown/kelvin_helmholtz/src/8b3059ee67b5db915c43f413d3d3e4e5fb3e3189/equations.py?at=default#cl-43

where the add_int_bc occurs. 

There may be some other tricks with the QG eqns, since I didn't see an obvious gauge BC in your script; can take a better look at it tomorrow.

--Ben

--
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/ECCB6097-D086-41DD-8412-43581C445CBB%40princeton.edu.
For more options, visit https://groups.google.com/d/optout.

[jsq...@pppl.gov]

Thanks for the quick reply. I realize now that my problem with the periodic domain was exactly what you say -- I needed a gauge condition. This came up on a previous formulation I had that used pressure, with the vorticity formulation it works fine without any extra condition (I had simply forgotten to check in the new script).

Thanks a lot!
Jono

On Monday, December 15, 2014 6:01:28 PM UTC+13, Ben Brown wrote:
> Jonathan,
>     Welcome!
>
>
> To answer your ps question: yes, you can do doubly periodic, but if you have a gauge condition you'll need to add an internal BC.
>
>
> There's an example here:
>
>
> https://bitbucket.org/bpbrown/kelvin_helmholtz
>
>
> See especially KH_periodic.py and
>
>
> https://bitbucket.org/bpbrown/kelvin_helmholtz/src/8b3059ee67b5db915c43f413d3d3e4e5fb3e3189/equations.py?at=default#cl-43
>
>
> where the add_int_bc occurs. 
>
>
>
> There may be some other tricks with the QG eqns, since I didn't see an obvious gauge BC in your script; can take a better look at it tomorrow.
>
>
> --Ben
>

[Daniel Lecoanet]

Hi Jono,

Sorry for the late response, but here's one way to add forcing in Dedalus.

Say you have a function f(t) which returns the forcing in grid space at a time t (you could pick f(t) to be white noise, or something more complicated).  Define a function as:

def forcing(solver):

    return f(solver.sim_time)

Then define a GeneralFunction:

forcing_func = GeneralFunction(domain,'g',forcing,args=[])

Then add a parameter to your problem called forcing_func, include the forcing in the equation, and point the parameter to the GeneralFunction:

problem = de.ParsedProblem( ... , param_names=['R','beta','forcing_func'])

problem.add_equation("dt(Z) - R*(dy(dy(Z)) + dx(dxZ))    = - dxP*(dy(Z)+beta) - dy(P)*dxZ + forcing_func")

problem.parameters['forcing_func'] = forcing_func

After you create your solver, you need to tell forcing_func that it's argument is solver:

forcing_func.args = [solver]

forcing_func.original_args = [solver]

This is a bit complicated, but we're working on simplifying the process...

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/befdff6a-d97b-47e2-87d8-0dd01a8163b2%40googlegroups.com.

[Jonathan Squire]

Thanks Daniel, I'll give this a go in a few days. 

One question, In the case that forcing_func is white noise, how will the solver know to scale the noise with sqrt(dt) rather than dt (as for deterministic forcing)?

Cheers,

Jono

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/BqTjYZzqHHw/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/CAJoYf%3Dhx8Lg_xwEG3D9uXkO_JNEqSw%2BY_Se3%3D9o5BF1NmQzPkQ%40mail.gmail.com.

[Daniel Lecoanet]

Good question.  You can make forcing be a function of dt, which you can access using solver.dt.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/3487E63F-1065-43E6-AECE-05E9D33F0034%40princeton.edu.

[Ben Brown]

Quick note:

    If you're using code post pull request #25 (generalized CFL), I think you may need to explicitly pass dt (or add it to the forcing method). I think solver no longer has a .dt attribute.

But I may be remembering wrongly...

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAJoYf%3DjbxkV5vWWaxL0AuDO%3DwtgPB8-UDfKL%3DuYsO5mQYd8KBg%40mail.gmail.com.

[Keaton Burns]

Yep I think that’s right.  I’ve added a quick docstring to GeneralFunction to make its usage a little clearer, but it’s really meant to be a light-weight wrapper, with any detailed logic being handled by the function underneath.

-Keaton

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

[joy.m...@gmail.com]

Just curious, but I notice from the code that solver.state.fields seems to contain the arrays containing the data. Can we use that as well in this formulation (i.e, return f(solver.state.fields)). I don't see any obvious reason why not!

Thanks,

Joy

[jsq...@clear.net.nz]

Thanks for all your help! That's exactly what I was looking for and everything seems to be working well now.

Best,
Jono

[Jack]

Hi Daniel, 

Has this process of adding time-dependent forcing been simplified? 

I try to do something similar. The syntax in this thread seems to be from an old version of Dedalus. 

 

Jack

> To unsubscribe from this group and stop receiving emails from it, send an email to dedalu...@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/ECCB6097-D086-41DD-8412-43581C445CBB%40princeton.edu.
>
> For more options, visit https://groups.google.com/d/optout.

--
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 dedalu...@googlegroups.com.

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