Using polar coordinates with UFL and C++

Question

Hello everybody,

I want to solve the following variational formulatiion using fenics:

$$\sum\limits_{i \in {h,s,b}} \int_{\Omega_i}\rho c_v \frac{\partial{\theta}}{\partial{t}}\,v\,dx - \int_{\Omega_i} \frac{\kappa_{r i}}{r} \frac{\partial \theta}{\partial r} \cdot v dx + \int_{\Omega_i} \kappa_{r i} \frac{\partial \theta}{\partial r} \frac{\partial v}{\partial r}+\frac{\kappa_{\phi i}}{r^2} \frac{\partial \theta}{\partial \phi} \frac{\partial v}{\partial \phi} \, dx$$

$$+ \int_{r=r_b}\alpha_b(\theta-\theta_{atm})\,v \, dS = 0$$

where $\Omega_{h,s,b}$ are some subdomains of a circle. My results with fenics are not like I expected them. While I was trying to look for mistakes, I wondered if I even use the right UFL-formulation? Is the following code the correct translation of the weak formulation into UFL-language, especially is this the right way to treat the
$\frac{1}{r}$?

V = FiniteElement("Lagrange", triangle, 1)

V_kappa = FiniteElement("DG", triangle, 1)

u = TrialFunction(V)
v = TestFunction(V)

dt = Constant(triangle)
u0 = Coefficient(V)
f = Coefficient(V)
rho = Coefficient(V_kappa)
cv = Coefficient(V_kappa)
kappa_r = Coefficient(V_kappa)
kappa_phi = Coefficient(V_kappa)
alpha = Constant(triangle)
u_atm = Constant(triangle);
x = triangle.x

Define bilinear and linear forms

a = rho * cv * (1/dt) * inner(u, v) * dx + alpha * u * v * ds(1) + (kappa_r * Dx(u,0) * Dx(v,0) + kappa_phi * 1/x[0] * Dx(u,1)* 1/x[0] * Dx(v,1)) * dx - kappa_r * 1/x[0] * Dx(u,0) * v * dx
L = rho * cv * (1/dt) * inner(u0, v) * dx + inner(f, v) * dx + alpha * (u_atm) * v * ds(1)

I am just trying to exclude possible error sources :)
Thanks!

Comments

It's the same approach I used for spherical and cylindrical coordinates in 1-D, where I just had to deal with r.

I would expect r=x[0], theta=x[1], and phi=x[2], with the spatial partial derivatives present in the forms:
Dx(x[1],x[2]) for example.

---

So it should be correct right? As a second approach I tried to manually use a coefficiention function
$$ pos = Coefficient(V), $$
used this function instead of 1/x[0] in the UFL-formulation and manually interpolated the expression given by $expression(x[0],x[1]) = 1/x[0]$ to pos, but the results differ a lot from the results given by the method of the original post.

Does anyone know how this is 1/x[0] is implemented? Because my first coordinate lies in the intervall $[0,R]$, but I dont get division by zero errors.

Thanks a lot!

---

To my understanding having x[i] is beneficial as the spatial coordinates are incorporated in the quadrature integration within the generated code. With a Coefficient(V), which is interpolated to x[i], the coefficient is integrated - introducing interpolation errors. I've tried it, it works, just not as well. I really only noticed a significant difference when I needed to find the divergence near r=0 in spherical coordinates (a great simplified problem for exploring this, find the div(r^2) near zero). Using x[i] performed much better in this case.

In my case I only solved from [eps, R], where I set eps to something reasonably small. As you pointed out, solving from [0, R] gave results, just nonsensical for the most part. I suspect some form of stabilization in the solving of the system with derivatives going to zero as r goes to zero .

Hopefully someone can chime in with more - I would love to have a nicer explanation.