Is there any efficient way to add many dirac deltas on a variational form?

1

Hi Fenics community,

I want to solve a non-linear version of the Poisson equation with point charges.

Becasue I am using the fenics NonLinearProblem class, I cannot use the PointSource function.

In here (Implement point source in the variational formulation), a single point source is included using a Dirac delta function.

I was wondering if is there any efficient way of do it with a huge number of point sources (on the order of 10^4)

Cheers,

Ximo

2

You can apply many point sources at once using this.

Furthermore, you can indeed use the NonlinearProblem class. Implement your own solver similar to here, but apply your point sources when your construct the residual (assuming they don’t constitute a component of the Jacobian).

3 Likes
3

I’ve used the approach that @nate suggested for large numbers of point forces corresponding to immersed boundaries before, and it’s decently efficient.

I’ll add that, if the point sources do contribute to the Jacobian, you can sometimes get away with a lumped-mass approximation of that contribution. For instance, if you have a mass-matrix-like point contribution

+\int C\delta(x-y)u(x)v(x)\,dx

to your variational form, where u and v are trial and test functions, then you can do something like

ps = PointSource(V,y,C)
lmVec = Function(V).vector()
ps.apply(lmVec)
as_backend_type(A).mat().setDiagonal(as_backend_type(lmVec).vec(),
                                     addv=PETSc.InsertMode.ADD)

where A is the LHS matrix from the rest of the variational form and PETSc is imported from petsc4py. This is not the exact Jacobian contribution, but can be good enough to prevent Newton iteration from diverging, and, if you converge with a consistent residual, the solution is unaffected by approximations in the Jacobian.

3 Likes
4

Hi Nate,

I also have the same problem with a large number of point sources in nonlinear problem. Could you please give a simple example about how to use PointSource class with nonlinear problem? Since as far as I remember, you need to do

A,b = assemble_system(a, L)

and

pointsources.apply(b)

which works for the linear problem?

Thanks a lot in advance.

Kind regards,
Alice