slightly complicated initial condition

[markadams]

I started with the Burgers equation demo but need an initial condition that is complex enough that I would much prefer not stuffing it into a string. https://www.firedrakeproject.org/interpolation.html#python-expression-classes looks fine to me but it has been removed. What is the best way to do this?

[dham]

Well we also don't support string expressions, so that's a happy coincidence ;).

If you've got an analytic expression for an initial condition then the preferred way to set the initial condition is to express that analytic experession in UFL and interpolate that. The key tools you need beyond the usual arithmetic and trig operations are SpatialCoordinate(mymesh) which will return a symbolic vector which you use as the coordinates in your expresion, and ufl.Conditional which you can use for any branches you need.

[markadams]

Thanks, I now need Atan2 (or some way to get a theta, and will need Bessel functions next. I see https://fenics.readthedocs.io/projects/ufl/en/latest/api-doc/ufl.html#ufl.classes.Atan2 but I get theta = Atan2(x[0],x[1]) NameError: name 'Atan2' is not defined

[rckirby]

UFL interpolation/projection — the things commented as “UFL equivalent” works.  It supports conditionals and basic transcendental functions as well.On Dec 19, 2022, at 6:25 AM, markadams ***@***.***> wrote: I started with the Burgers equation demo but need an initial condition that is complex enough that I would much prefer not stuffing it into a string. https://www.firedrakeproject.org/interpolation.html#python-expression-classes looks fine to me but it has been removed. What is the best way to do this? —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[rckirby]

`atan_2`, `bessel_I`, `bessel_K`, `bessel_J`, and `bessel_Y` all enter the namespace after `from firedrake import *` Hopefully these will DTRT for you

…

On Mon, Dec 19, 2022 at 8:40 AM markadams ***@***.***> wrote: Thanks, I now need Atan2 (or some way to get a theta, and will need Bessel functions next. I see https://fenics.readthedocs.io/projects/ufl/en/latest/api-doc/ufl.html#ufl.classes.Atan2 but I get theta = Atan2(x[0],x[1]) NameError: name 'Atan2' is not defined — Reply to this email directly, view it on GitHub <#2695 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4XUL454WCNNSEW5VJWRPLWOBXUNANCNFSM6AAAAAATDJS5ZQ> . You are receiving this because you commented.Message ID: ***@***.*** com>

[markadams]

Great, reproduced the paper's IC, Thanks.
And let me know if you have any 'best practices' suggestions.
(this is an MHD equilibrium)

x = SpatialCoordinate(mesh)
r = sqrt(dot(x,x))
theta = atan_2(x[0],x[1])
cos_theta = cos(theta)
outer_f =  (1/r - r)*cos_theta
inner_f = (2/(Constant(k)*bessel_J(0,Constant(k))))*bessel_J(1,Constant(k)*r)*cos_theta
f = conditional(gt(r, 1.0), outer_f, inner_f)
ic = project(f, Psi)

[wence-]

The only best practice with complicated (non-polynomial) expressions like this that you're missing is likely manual specification of the quadrature degree (so that we don't estimate some ludicrous quadrature for the integration). You might also want to consider interpolateing your initial condition rather than projecting. I doubt it matter in this case, but if you have certain analytic properties of an initial condition (e.g. divergence-free-ness) then the L2-projection is not guaranteed to preserve those whereas interpolation will (up to the accuracy of the moments in you finite element's functionals).