Point source of electric current in radially symmetric cylindrical coordinates

Hello, Netgen/NGSolve community,

I am trying to model the distribution of electric potential around a point source of DC electric current in radially symmetric cylindrical coordinates.

The problem I am trying to solve is given by the following equation:
[tex]\nabla \cdot (\sigma \nabla u) = - \nabla \cdot j[/tex]
with the source term
[tex]\nabla \cdot j = I \delta (x-x_s) \delta (y-y_s) \delta (z-z_s)[/tex]
where sigma is the distribution of conductivity in the domain, u is the electrostatic potential, j is the electric current density, I is a current injected at the point source located at (x_s, y_s, z_s) and delta is the Dirac’s delta.

In radially symmetric cylindrical coordinates (r, z), with the source located at (0, 0) this equation can be rewritten as:
[tex]\frac{1}{ r} \frac{ \partial}{ \partial r} [ \sigma \frac{\partial u}{\partial r}] + \frac{ \partial}{ \partial z} [ \sigma \frac{\partial u}{\partial z}] = \frac {I \delta (r) \delta (z)}{2 \Pi r}[/tex]

Originally I was using GMSH and FEniCS, but over time I realised that switching to Netgen/NGSolve might be beneficial since model geometry I am dealing with is easier to describe in Netgen than in GMSH. In addition, this will allow me to skip mesh conversion from GMSH format to FEniCS format which is surprisingly time-consuming.

Bellow, I supplied the code I managed to create. In the final implementation the domain around the current source will be much more complex but for the sake of simplicity of an example and to allow comparison with the analytical solution let’s assume that the current source is located in a homogenous medium. At the end of the code I’ve added a comparison with an analytical solution:
[tex]U (r_I) = \frac {I} {4 \Pi r_I \sigma } [/tex]
where r_I is a distance between the source point and the point at which the potential value is measured.

The numerical and analytical results should match (taking aside drop of potential caused by the finite size of the domain), but unfortunately, they are pretty far apart in case of my Netgen/NGSolve implementation while in FEniCS similar code gives me correct results. Since I have no experience with Netgen/NGSolve I probably messed something up. I guess that the problem is associated with the source term. In FEniCS I used build-in function PointSource() to add point source in one point of the domain. To my knowledge similar function is not present in Netgen/NGSolve, therefore I added point source using a coefficient function with Dirac’s delta approximation.

I would be very grateful if someone will help track the source of error. In addition, I wonder if it is possible to assign maxh value to the point or project it as a coefficient function on the entire mesh. In all examples I was looking at the maxh value was assigned to lines. In my implementation, I quite crudely added two short lines on both sides of the source location and assing maxh value to them.

EDIT: Attachments are in zip format because I was unable to add attachments as individual images.

Hello,

welcome to NGSolve !

I put two versions for the point-sources into the attached file:

1.
fixing the approximation by the Gaussian peak: you need the radius-factor also for the right hand side.
In the file we are checking that the integral approximates 1.

2.
An AddPointSource Python-function.
It searches the element, evaluates shape functions, and adds the contributions to the linear-form.
It must be called after f.Assemble()


As you give the maxh to the geometry edges, you can also assign a maxh to the geometry points

Best,
Joachim

Thanks a lot, both versions work perfectly.

we have now a built in version of point sources (available in the current nightly):

best, Joachim

Dear Joachim,

If I am solving a vector problem using HDG. How can I include pointwise sources using the function ?

The V space is defined as:

And the source is characterized as a vector field with intensity I, given by:
[tex]
\begin{equation}
\mathbf{I} = \left(
a + b \mathbf{i},
c + d \mathbf{i}
\right)^T
\end{equation}
[/tex]

[tex]
\begin{equation}
\mathbf{f} = \mathbf{I}\delta(x-x_0)
\end{equation}
[/tex]


Ultimately, I'd like to introduce several pointwise sources
[tex]
\begin{equation}
\mathbf{f}(\mathbf{x}) = \sum_{j = 1}^{n} \mathbf{I}_j \delta(\mathbf{x} - \mathbf{x}_j).
\end{equation}
[/tex]

But, CompoundFiniteElement has no attribute 'CalcShape'.

Is there a way to use the function to generate the pointwise sources?


Best regards,

Alan.