Cylinder symmetric Poisson solution

[Andreas Müsing]

Dear All,

I am solving a cylinder symmetric Poisson problem in deal.ii. 

My program is already working properly, but I found the solution by trial&error, and my hope is to get the final understanding of what I did via the user group.

My question is closely related to  this stack exchange question  , wherethe math is done for an equation independent on the z coordinate, but depending on theta.

The Poisson equation (strong form) in cylinder coordinates with cylinder symmetry is

1r∂∂r(r∂∂rϕ(r,z))+∂2dz2ϕ(r,z)=−ρ(r,z)

In order to avoid the sinularity at r = 0, Prof. Bangerth suggested in the SE link to multiply with the proper area element 2 p r dr dz in order to get the weak form\frac{1}{r} \frac{∂}{∂r}(r \frac{∂}{∂r} \phi(r,z)) + \frac{∂^2}{dz^2}\phi(r,z) = - \rho(r,z)

∫0R∫−∞∞φ(r,z){∂∂r(r∂∂rϕ(r,z))+r∂2∂z2ϕ(r,z)}drdz= -∫0R∫−∞∞φ(r,z)ρ(r,z)rdrdz

After integration by parts, I get the following equation (where I am not so sure anymore if it is correct):

∫0R∫−∞∞r∂φ∂r∂ϕ∂r+r∂φ∂z∂ϕdzdrdz=∫0R∫−∞∞ϕρrdrdz\int_0^R \int_{-\infinity}^\infinity r \frac{∂\varphi}{∂r} \frac{∂\phi}{∂r} + r \frac{∂\varphi}{∂z}\frac{∂\phi}{dz} dr dz = \int_0^R \int_{-\infinity}^\infinity \phi \rho r dr dz

My questions now are

a) Is the weak form above correct for the cylinder symmetric case?

b) How do I implement this in deal.ii? 

c) When I look at the first term in the PDE's strong form, I can carry out the first term r derivative by applying the derivative product rule, which results in terms having a first and a second derivative in r, respectively. Therefore, the solution should be in any case a linear combination of first derivative and second derivative values.

I used tutorial step-63 with success as a starting point, where the advection-diffusion equation has both a second derivative (diffusion) term, and a first derivative (=advection) term. Further, I used an analytic known solution of my problem to compare the results, and I get reasonable solutions for e = 1 and b = 2. But why beta = 2 and not beta = 1?

Another way of looking at the problem: I found just by experimenting, that the following to deal.ii implementations give the same (reasonable) results for my problem:

//--- first implementation ----

double advection_direction[0] = 2.0, double advection_direction[1] = 0;

cell_matrix(i, j) +=
                fe_values.JxW(q_point)* (

                (radius * fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)) +
                (fe_values.shape_value(i, q_point) * (advection_direction * fe_values.shape_grad(j, q_point) );

cell_rhs(i) += radius * radius * some_value;

//---- second implementation, but identical results ---

double advection_direction[0] = 1.0, double advection_direction[1] = 0;

cell_matrix(i, j) +=
                fe_values.JxW(q_point)* (

                (fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)) +
                (1.0 / radius * fe_values.shape_value(i, q_point) * (advection_factor * fe_values.shape_grad(j, q_point) );

cell_rhs(i) += some_value;

The difference between both versions seems to be a factor r² in the weak form. But I don't understand why the advection term needs a different value ( 1 versus 2) for both versions. My impression is that the factor of 2 is a result of the integration by parts, when including a factor r² (its derivative is 2r), but I don't see the actual connection between my experimental results and the math above.

Any help or hint into the right direction is appreciated

Regards
Andreas

[Andreas Müsing]

P.S: Please find in the attachment a screenshot of the equations, they were unfortunately removed/scrambled in my original post.

[Wolfgang Bangerth]

On 4/23/24 07:59, Andreas Müsing wrote:
>
> a) Is the weak form above correct for the cylinder symmetric case?

Yes, this looks correct. It also correctly leads to a symmetric bilinear form,
which is reassuring for the laplace equation.

> b) How do I implement this in deal.ii?

for (cell=...)
for (q=...)
for (i=...)
for (j=...)
local_matrix(i,j) += fe_values.shape_grad(i,q) *
fe_values.shape_grad(j,q) *
fe_values.quadrature_point(q)[0] * // r
fe_values.JxW(q);

> c) When I look at the first term in the PDE's strong form, I can carry out the
> first term r derivative by applying the derivative product rule, which results
> in terms having a first and a second derivative in r, respectively. Therefore,
> the solution should be in any case a linear combination of first derivative
> and second derivative values.

Can you elaborate? It's true that if you multiply everything out, the equation
is of the form
d/dr ... + d^2/dr^2 ... = ...
but I don't understand how you infer something about the solution then?

> I used tutorial step-63 with success as a starting point, where the
> advection-diffusion equation has both a second derivative (diffusion) term,
> and a first derivative (=advection) term. Further, I used an analytic known
> solution of my problem to compare the results, and I get reasonable solutions
> for e = 1 and b = 2. But why beta = 2 and not beta = 1?

What is the form of the advection equation you are trying to solve, in
cylindrical coordinates?

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Andreas Müsing]

Wolfgang Bangerth schrieb am Dienstag, 23. April 2024 um 16:45:36 UTC+2:

On 4/23/24 07:59, Andreas Müsing wrote:

> b) How do I implement this in deal.ii?

for (cell=...)
for (q=...)
for (i=...)
for (j=...)
local_matrix(i,j) += fe_values.shape_grad(i,q) *
fe_values.shape_grad(j,q) *
fe_values.quadrature_point(q)[0] * // r
fe_values.JxW(q);

this was the first approach that I tried, but the results were not what I expected. My answer to your reply on c) will clarify it.

 

> c) When I look at the first term in the PDE's strong form, I can carry out the
> first term r derivative by applying the derivative product rule, which results
> in terms having a first and a second derivative in r, respectively. Therefore,
> the solution should be in any case a linear combination of first derivative
> and second derivative values.

Can you elaborate? It's true that if you multiply everything out, the equation
is of the form
d/dr ... + d^2/dr^2 ... = ...
but I don't understand how you infer something about the solution then?

I have a magnetostatic model, where I calculate the magnetic vector potential A(r,z) in cylinder coordinates, assuming a cylinder symmetric current density J(r,z).

For a delta peak excitation in J(r,z), the magnetic field B = curl A is known analytically on the symmetry axis (B(r=0, z)). This delta excitation (in FEM parlance: the rhs vector has a single non-zero component) is a current circle filament in 3D space.

For the same current circle filament, I created a separate verification solution via the Biot-Savart law over the whole domain, to get an alternative solution B(r,z).

Your implementation suggestion (which I tried before), introduces a factor r in the matrix. This by itself did not lead to the expected spatial solution of A(r,z) and B(r,z), when I compared with B_z(r=0) or the Biot-Savart solution.

Then my reasoning was, that the involved integrals contain the expression ∂/∂r (r ∂A / ∂r) =  ∂A/∂r + r ∂^2/∂r^2 A , and I assumed that the Laplace equation in cylinder coordinates requires an additional first derivative term.

 

> I used tutorial step-63 with success as a starting point, where the
> advection-diffusion equation has both a second derivative (diffusion) term,
> and a first derivative (=advection) term. Further, I used an analytic known
> solution of my problem to compare the results, and I get reasonable solutions
> for e = 1 and b = 2. But why beta = 2 and not beta = 1?

What is the form of the advection equation you are trying to solve, in
cylindrical coordinates?

There was a misunderstanding: I am still solving a Laplace problem, but the transition from Cartesian coordinates to cylinder coordinates added a term to the equation which can be interpreted as an advection term. The structure of the PDE in step-63 is therefore exactly what I need. I played with the parameters epsilon and beta, and for the mentioned epsilon=1.0 and beta = 2.0, the solution matched exactly my analytic results.

I am happy with my implementation except of the important question that bothers me: why beta = 2? Maybe I did somewhere else a silly mistake, but most probably  I just don't have yet enough insight into FEM and its math details.

Thanks for your help!

[Andreas Müsing]

Hello  Wolfgang

you last, I reconsidered my original solution and indeed found a silly mistake when I calculated the curl of the potential values. Of course then, the implementation that you suggested was correct and simple enough.

Probably, when I added the additional first derivative term, I solved the equation in r * A (or radius * potential), which was fixing the bug in my evaluation of the B field.

Sorry for my email flood, and thanks again for your help, hopefully this topic will be useful for someone else in the future.

Andreas