How to Implement a Point Source with the Weak Form

Chien Liu August 24, 2015
Share this on Facebook Share this on Twitter Share this on Google+ Share this on LinkedIn

Today we continue our discussion on the weak formulation by looking at how to implement a point source with the weak form. A point source is a useful tool for idealizing the situation where a source is concentrated in a very small region of the modeling domain. We will find that it is very convenient to set up such a point source using the weak form.

The Mathematics of a Point Source

Consider a one-dimensional domain on the x-axis with a source localized around x = 0. We can plot the strength of the source as a function of x and it may look like this:

Plot of a one-dimensional domain localized around a source.

Here, we have assumed that the strength has a constant value of 1/w within the interval [-w/2, w/2] and is zero everywhere else. This gives a rectangular shape of width w and height 1/w, as shown in the figure above. The function is often called a rectangular, top-hat, or sometimes, a disc function. The total strength of the source is given by the area of the rectangle, which is unity.

For linear systems, if we only care about what happens far away from the source where \left| x \right| \gg w, then the actual shape of the source strength does not matter much, as long as the area beneath that shape is the same. Furthermore, we are free to make w progressively smaller and smaller: the width of the rectangle decreases while its height increases in such a way that the total area remains the same, as shown in the graph below.

Schematic showing change in the localized source.

The localized source represented by the blue curve is progressively made thinner and taller (the orange and green curves), while maintaining the integrated strength of unity.

Eventually, we arrive at a rectangle that is infinitesimally thin and infinitely tall, but still has a well defined area of unity. This leads us to the so-called delta function \delta(x) and, correspondingly, the localized source now becomes an idealized point source of unit strength.

The delta function has some convenient properties. Its value is zero everywhere except at the origin:

\delta(x)= \left\{ \begin{array}{ll}
\infty \mbox{ for } x=0\\
0 \mbox{ elsewhere}
\end{array} \right.

Integrating the product of a delta function and another function just extracts the value of the latter function at the origin:

\int^\epsilon_{-\epsilon} \delta(x) f(x) dx=f(0) \mbox{ for all } \epsilon > 0

A point source at a general position x=a can be obtained by a simple coordinate shift of the delta function \delta(x-a). We have

\delta(x-a)= \left\{ \begin{array}{ll}
\infty \mbox{ for } x=a\\
0 \mbox{ elsewhere}
\end{array} \right.

and

\int^{a+\epsilon}_{a-\epsilon} \delta(x-a) f(x) dx=f(a) \mbox{ for all } \epsilon > 0

It is also easy to generalize the delta function and the corresponding point source to higher dimensions. For example, in 2D, we have

\delta(x-a,y-b)= \delta(x-a)\delta(y-b) = \left\{ \begin{array}{ll}
\infty \mbox{ for } x=a, y=b\\
0 \mbox{ elsewhere}
\end{array} \right.

and

(1)

\iint\limits_{\Omega} \delta(x-a,y-b) f(x,y) dx dy=f(a,b) \mbox{ for all } \Omega \ni (a,b)

Implementing a Point Source Using the Weak Form

This tutorial solves the Poisson equation on a unit disc with a point source at the origin. The equation reads

(2)

-\left({\partial^2 \over \partial x^2}+{\partial^2 \over \partial y^2} \right)u(x,y) = \delta(x,y)

where u is the dependent field variable to be solved.

At first sight it may not be obvious how to discretize this equation to be solved numerically. What value do we put at the origin for the source term on the right-hand side? The value of the delta function is infinite there, but computers don’t like infinities!

Here, we will see that the weak formulation comes in handy. Recall that in this introductory blog post on the weak form, we multiply the differential equation to be solved by a test function and integrate over the entire domain (See Eq.(4) in that post). We can follow the same procedure here to solve Eq. (2). After multiplying by a test function \tilde{u}(x,y) and integrating over the unit disc domain, the right-hand side of Eq. (2) simply becomes

(3)

\iint\limits_{\Omega} \delta(x,y) \tilde{u}(x,y) dx dy=\tilde{u}(0,0)

by using the integration property of the delta function given in Eq. (1). This gives us something very easy to implement in COMSOL Multiphysics.

Start with a new 2D model with the Weak Form PDE physics interface and a Stationary study. Draw a unit circle centered at (0,0) and draw a point there as well. Set the Weak Expressions field under the default Weak Form PDE 1 feature to -test(ux)*ux-test(uy)*uy. This takes care of the left-hand side of Eq. (2) in exactly the same fashion as for the 1D case discussed in this previous post.

Now, for the point source on the right-hand side, \tilde{u}(0,0), we simply add a point Weak Contribution node and select the point at the origin. For the Weak expression, we enter test(u). It’s that simple for the point source!

Screenshot showing how to add a Weak Contribution node in COMSOL Multiphysics.
It may be worth noting that by entering test(u), we set the strength of the point source to unity. For any other source magnitude, simply multiply by a factor. For example, the expression 2*test(u) gives a point source of strength 2.

After finishing the set-up with a Dirichlet boundary condition at the perimeter of the circle, we can solve the model and observe the same solution as seen in the point source tutorial mentioned above:

Schematic of the point source tutorial results.

Also as seen in the tutorial, the numerical solution (blue curve) matches the analytical solution (green curve) very well, except near the original where a singularity occurs:

Graph comparing the numerical and analytical solutions.

As mentioned earlier, the point source provides a convenient idealization of a localized source in situations where we only care about the solution far away from the source. We illustrate this point with the following graph, where we have added three more curves to the graph above. These three curves are numerical solutions to the same Poisson equation in the same unit disc domain, but with various sizes of top-hat, or disc, shaped sources replacing the point source. The integrated strength of each top-hat source is calibrated to unity by setting its height to one over its area, in the same fashion as in the 1D case shown in the image above. As we see clearly from the figure below, all solutions are indistinguishable from one another far away from the sources. (In this example for x \gg 10 \, mm.)

Plot comparing different point source solutions and disc source measurements.

Conclusion

Here, we have demonstrated the ease of creating point sources using the weak form. The numerical difficulty in the representation of the delta function is circumvented with a simple integration. In upcoming posts we will look at discontinuities and boundary conditions. Stay tuned!


Categories


Comments

  1. Vishwas Nesarikar August 25, 2015   12:17 am

    Very nice explanation. Thanks

  2. Perry Gao November 11, 2015   4:04 am

    Dear Chien,

    Thanks for your great post!
    I recently met some problem in implementing ‘identity pairs’ to realizing discontinuous boundary conditions in weak form pde.
    Is this going to appear in your next blog?
    Thanks!

    Perry

  3. Carolina Montoya-Pachongo November 28, 2016   12:47 pm

    Dear Chien, wonderful posts about the weak form of PDEs. I can see you did not publish posts anymore. Is someone else in charge of the following topics: discontinuities and boundary conditions? Many thanks.

  4. Chien Liu November 28, 2016   1:47 pm

    Dear Carolina, yes we do hope to continue the blog series in the future. Stay tuned! Chien

Loading Comments...

Categories


Tags