# How to implement point source or volume source in finite element implementations

I'm trying to do a simple implementation to study the advection-diffusion-reaction dynamics in a straight pipe.

I have points positioned along the length of the pipe (blue dots in the image above).

I came across a blog (link) that explains how a point source is implemented. However, I am not really sure how to adapt it in this case.

The equation that captures the concentration dynamics is the following

The R in the last term of RHS corresponds to sink/source and the units is generally defined in mol/m^3 s. I want to define R = kc where the units of k is per sec. But I want to define R only at the points colored blue (i.e only as a point source/sink and not as a line source or sink).

I would like to ask for suggestions on how this implementation can be done.

EDIT: I would like to compare the FEM result obtained for 1D straight pipe simulation of convection-diffusion dynamics with point sinks located along the length of the pipe (problem posted above) against the analytical solution.

Suggestions on the analytical expression against which the numerical result can be compared will be really helpful.

• If you want to add them as point sources in your nodes I think that you can add them to the right have side vector. Jun 2 at 13:33
• @nicoguaro From what I understand, by default software like COMSOL allows the implementation of a line source since the source/sink term R is associated with area. Could you please elaborate? I'm a bit confused about how to add a point source/sink here. Jun 2 at 14:14
• Does the documentation of deal.II's VectorTools::create_point_source_rhs() function help? dealii.org/developer/doxygen/deal.II/… Jun 2 at 14:29
• You could have something like $R = \delta(x - x_i)$ to represent a point source. Jun 2 at 20:04
• I think that should be another question. Jun 22 at 3:01

Seen that you also have an account over at mathematica stackexchange I am going to show an implementation using Mathematica. This is certainly not the only way you can do it but hopefully gets you started.

We start by creating a 1D mesh

Needs["NDSolveFEM"]
region = Line[{{0}, {1}}];
includePoints = {{1/3}, {2/3}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints,
"MaxCellMeasure" -> 0.0008]


The mesh has points included at 1/3 and 2/3. This is where we will position the point sources later. I also refined the mesh, because I am going to add a convection term which will make this PDE convection dominant (large Peclet number).

We set up the dependent variable, the mass concentration $$c$$, time variable $$t$$ and the spatial independent variable $$x$$. Also some material properties are defined.

vars = {c[t, x], t, {x}};
pars = <|"DiffusionCoefficient" -> 1,
"MassConvectionVelocity" -> {100}|>;


For the point sources the idea is to use a regularized delta function. The Dirac delta function poses a problem in numerical simulations with FEM implementation that is not prepared to take it as an input. This is because the Dirac delta function is singular at the source location $$X_s$$ and needs special treatment that differs from how FEM treats other terms. This can be avoided with an approximation to the Dirac delta function. The process of approximating the Dirac delta function is called regularization.

RegularizedDeltaPoint[g_, X_List, Xs_List] :=
Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]


More information about this can be found in the Heat Transfer tutorial. The process for mass transport is the same. So if you also want to volume sources you can proceed in the same way as mentioned in that tutorial.

Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;


The regularization function needs a parameter gamma, which is based on the minimal element diameter. Let's construct the point source at the first included point and visualize it.

temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x},
includePoints[[1]]];
Plot[temp, {x, 0, 1}]


You see that this is an idealized approximation to the delta distribution. Note that

Integrate[temp, {x, 0, 1}]
1.


The Comsol approach uses a weak formulation. That is a different approach. I am trying to show that (point) sources can be set up in different ways.

Qp = 10;
pars["MassSource"] =
RegularizedDeltaPoint[Subscript[gamma, reg], {x},
includePoints[[1]]]*Qp +
RegularizedDeltaPoint[Subscript[gamma, reg], {x},
includePoints[[2]]]*If[t > 1/2*10^-3, 5*Qp, 0];


Qp is the point source strength. The unit would be [$$mol/(m^3 s)$$], because the regularization function is dimension less the source has the same units as the overall PDE components, see here. The first point source is active all the time while the second only starts at a specific time but is stronger.

We set up the PDE and solve it.

pde = {MassTransportPDEComponent[vars, pars] == 0, c[0, x] == 0};
tEnd = 10^-3;
cfun = NDSolveValue[pde, c, {t, 0, tEnd}, {x} \[Element] mesh];


Visualize / explore the result.

Manipulate[
Plot[cfun[t, x], {x} \[Element] region,
PlotRange -> {{0, 1}, {0, 1/2}}], {t, 0, tEnd}]


More information can be found in the Mass Transport tutorial and the MassTransportPDEComponent.

I hope this is useful to you.

Update for older version:

The code I posted above works for version 12.3. Since the OP has 12.0 the following is the same code, suitable for version 12.0. In essence the MassTransportPDEComponent is not available in 12.0. But that's not a problem as one can write the PDE manually. Here is the complete code:

This first part is the same as before:

RegularizedDeltaPoint[g_, X_List, Xs_List] :=
Piecewise[{{Times @@
Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}];

Needs["NDSolveFEM"]
region = Line[{{0}, {1}}];
includePoints = {{1/3}, {2/3}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints,
"MaxCellMeasure" -> 0.0008];

Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;


This part is where things change. I added a parameters section and added a manual form of the PDE.

parameters = {kappa -> {{1}}, v1 -> 100,
gamma -> Subscript[gamma, reg], Qp -> 10};
pde = {Derivative[1, 0][c][t, x] +
Inactive[Div][(-kappa) .
Inactive[Grad][c[t, x], {x}], {x}] + {v1} .
Inactive[Grad][c[t, x], {x}] -
Qp*RegularizedDeltaPoint[gamma, {x}, {1/3}] -
If[t > 1/2000, 5*Qp, 0]*
RegularizedDeltaPoint[gamma, {x}, {2/3}] == 0,
c[0, x] == 0} /. parameters;

tEnd = 10^-3;
cfun = NDSolveValue[pde, c, {t, 0, tEnd}, {x} \[Element] mesh];


This should give you the same result in version 12.0.

Update: Cloud version

Another alternative is to use the wolfram cloud. Here is a cloud version of the notebook.

• Hi, I've copied the lines of code posted by you to try this out. Please find the notebook here. I get the following error :NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations.. Could you please help me in resolving this? Jun 18 at 13:30
• @Natasha, what version do you have? Jun 18 at 13:36
• I'm using version 12. Jun 18 at 13:39
• Depends a bit on how you want things to be formatted, you could start with this: data = Table[ Flatten[cfun[t, #] & /@ includePoints], {t, 0, tEnd, tEnd/10}]; Export["~/data.txt", data, "CSV"] Jun 21 at 4:55
• @Natasha, that's correct. I fixed it in the text. Jun 22 at 4:54

Contrary to what @user21's answer, I don't think that you need to do anything special for point loads. Let's see why.

A point load can be represented as a Dirac delta "function". So, in your case it would be something like

$$R = \rho \delta(x - x_i)\, ,$$

where $$\rho$$ is the intensity of the source and $$x_i$$ is the position.

If we use a weighted residuals for the finite element method you end up with the following integral for "body" sources

$$b = \int_\Omega R(x) w(x)\, ,$$

being $$w(x)$$ the weight function. If $$R(x) = \rho \delta(x - x_i)$$,

$$b = \int_\Omega \rho \delta(x - x_i) w(x) = \rho w(x_i)\, .$$

Thus, it suffices to evaluate your weight function in the point $$x_i$$. For "common" finite elements, we use interpolating functions. This implies the following

$$b = \begin{cases} \rho\quad \text{for nodes}\, ,\\ \rho \phi_j(r_i)\quad \text{otherwise}\, , \end{cases}$$

where $$r_i$$ represents the position of the point in the reference element.

Notice the following:

• For points that coincide with nodes (your case) the source contribute to exactly one entry in the $$b$$ vector. But for points inside elements it contributes to all the nodes related to that element, one for each function $$\phi_j$$.

• For linear elements $$\phi_j(r_i) > 0$$, but that's not the case for higher-order elements.

Regarding the question if one should or should not use point sources and its physical interpretation... It is a different matter, more related to the physics of the problem, I would argue. Respecting the units that you mention at the end of your question, I think that you don't provide enough information to assess it.

• As I said, there are many ways to do it. Jun 4 at 9:51
• @user21, yes. But you also said this: "The Dirac delta function poses a problem in numerical simulations as it can not be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location Xs. Hence, an approximation to the Dirac delta function is needed" Jun 4 at 16:19
• I see what you are saying, but that has the disadvantage that for a general purpose FEM code you'd need to introduce some spacial case code, no? - my approach works for other FEM codes that do not have this specifically build in. I have reworded my answer a bit to reflect that. Jun 7 at 4:46
• @nicoguaro Thanks a lot for the detailed explanation. I tried setting up the system described in my original post i.e. add point source at the blue points. As shown in the example demonstrated here I defined Q*test(u) in the weak form PDE. But the solution failed to converge. I am not sure if this is due to the singularity of the delta function at the source point. Jun 7 at 16:03
• @Natasha, I could imagine that this does not work in comsol, as when you specify the expression with the test function Q*test(u) comsol will still integrate over it, and, if I understand nicoguaro's answer correctly, you do not want this. Again, I do not have access to comsol and this is a shot in the blue. You could try to impose this as a DirichletCondition where you use comsol's h(x) u = r(x) and have h(x)=1 and r(x)=if(x==POS), Q, 0). maybe that works. Or, as an alternative, my approach will work. Specify the regularization function and then call it in the load vector. Jun 15 at 3:54