1
$\begingroup$

In my earlier question Interpolating the gradient of a cylindrically symmetric potential field expected obey the Laplace equation, especially near/across r=0 I give a simple example of a Jacobi relaxation script - the Laplace equation is solved for a cylindrically symmetric electrostatic lens system. The symmetry reduces the problem to 2 dimensions.

I am reading and learning about FEM and rather than write my own I would like to start with an existing package to solve the same cylindrical lens problem (ideally in 2D for speed) with all Dirichlet boundary conditions (i.e. I define the potential on all surfaces and boundaries).

I have seen PySfe, pyGIMLi, Fenics and Firedrake and scikit-fem Most if not all of these have good examples for Laplace in 2D but these are cartesian examples.

What I need is

  1. A recommendation for a package that is easy to install (I'm not a developer, I can spell "pip" and "conda" (barely) but not much more) and will have some way to generate solutions I can compare directly with my Jacobi relaxation work.
  2. Advice how to modify the (likely present) Laplace example for the problem as 2D cylindrical, i.e. $\phi(r, z)$ and $\theta$-independent.

Ideally I'll set up the problem in an IDE and run from a command line. I'm using macOS with an Anaconda installation.


From the linked question

Laplace equation Jacobi relaxation for cylindrical lens

$\endgroup$
2
  • $\begingroup$ Please provide the actual equations that you would like to solve and the necessary boundary conditions instead of a reference to your previous experimentation with this problem. $\endgroup$ Commented Feb 8 at 15:11
  • $\begingroup$ @Konstantinos if I'm not mistaken (it's 2:30 AM here) Laplace's equation expressed in cylindrical coordinates with no $\theta$ dependence will be $$\frac{\partial^2 \phi}{\partial r^2} + \frac{1}{r} \frac{\partial \phi}{\partial r} + \frac{\partial^2 \phi}{dz^2} = 0.$$ $\endgroup$
    – uhoh
    Commented Feb 8 at 18:34

1 Answer 1

2
$\begingroup$

You did not provide all the details so I made some assumptions. First, assuming symmetry in $\theta$, the Laplace equation in cylindrical coordinates reads $$\partial_r^2 u + \frac{1}{r} \partial_r u + \partial^2_z u = 0.$$ Let $\Omega = (0, 1)^2$ where the variables are $(r,z)$. The symmetry boundary is $\Gamma_0 = \{ (r, z) : r = 0 \}$ and the remaining part is $\Gamma_1 = \partial \Omega \setminus \Gamma_0$. We consider the boundary conditions $$\partial_n u|_{\Gamma_0} = 0, \quad u|_{\Gamma_1} = g$$ for given boundary data $g$.

For testing, an analytical solution is given by the Bessel function of the first kind $J_0$, see a question on another forum. In particular, $$J_0(r)(\exp(z) + \exp(-z))$$ satisfies the above Laplace equation when its restriction onto $\Gamma_1$ is used as $g$.

We define the function spaces $$V_z = \{ w \in H^1(\Omega) : w|_{\Gamma_1} = z \}$$ for a given $z$. The weak formulation reads: find $u \in V_g$ such that $$\int_\Omega -\partial_r u \partial_r v + r^{-1} \partial_r u v - \partial_z u \partial_z v = 0 \quad \forall v \in V_0.$$

After pip install scikit-fem==9.0.1 you can implement the weak formulation as follows:

from skfem import *
from skfem.helpers import *

@BilinearForm
def laplace(u, v, w):
    r, z = w.x
    ur, uz = u.grad
    vr, vz = v.grad
    return - ur * vr + 1 / r * ur * v - uz * vz

Then you can solve the problem with the boundary data from the Bessel function:

from scipy.special import jv

m = MeshTri().refined(4)
basis = Basis(m, ElementTriP1())
A = laplace.assemble(basis)

y = basis.project(lambda x: jv(0, x[0]) * (np.exp(x[1]) + np.exp(-x[1])))
x = solve(*condense(A, 0 * y, x=y, D=basis.get_dofs({'top', 'right', 'bottom'})))

basis.plot(x, shading='gouraud', colorbar=True).show()

If you have matplotlib installed, running the above code will result in the following figure:

Solution to the cylindrical Laplace equation.

You can also plot the difference x - y to see that the error is small:

Pointwise error.

I was a bit suspicious of the larger error at $r=0$ so I tried different mesh refinements to see that it converges as expected. I concluded that the error is probably due to the term with $1/r$ which warrants for more integration points or more elements near $r=0$.

(Note: I am the maintainer of scikit-fem.)

$\endgroup$
2
  • $\begingroup$ This is exactly what I need to get started. I spent yesterday at the library reading books and viewing online material but found myself saying "OK I believe you, very rigorous, but how do I do it?" A quick installation question; yesterday I'd already used pip install scikit-fem and pip install meshio (I have matplotlib from my anaconda installation) and while your example runs fine, skfem.__version__ yields 8.1.0 If I use your command for 9.0.1 and add ` -U` should that bring me up to date? $\endgroup$
    – uhoh
    Commented Feb 9 at 2:20
  • 1
    $\begingroup$ Yes, that should work. $\endgroup$
    – knl
    Commented Feb 9 at 4:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.