PointSource in dolfinx

dolfinx
1

Hello everyone,
I work with inhomogeneous Helmholtz equation and i wonder if the PointSource that was available in dolfin will be implemented in dolfinx.
Right now i am approximating it with the following steps:

  1. delta = interpolate narrow normal distribution
  2. integrate delta over volume
  3. defining PointSource = delta/integral(delta)

Do you agree that it is a correct way to do that?
How about the implementation of Pointsource in dolfinx? It would be way easier if it was implemented at the vector assembly level.

Thank you

2

That would be one way to approximate a point source. Even a ā€œtrueā€ point source will typically result in low-order accuracy anyway (due to the low regularity of the exact solution), so a Gaussian with O(h) standard deviation that is explicitly normalized using the same quadrature rule that it will be integrated with in variational forms is probably reasonable. (Using the same quadrature rule is important, especially if the standard deviation is small relative to the mesh size, since there could be significant quadrature error.)

Re-implementing PointSource for DOLFINx is still an open issue on GitHub, but has been open since 2018, so I donā€™t think itā€™s a very high priority, and I wouldnā€™t expect it to be resolved soon.

1 Like
3

I am planning to reimplement it next week. It should be quite straightforward now.

3 Likes
4

Amazing!
As always, thank you very much.
With that, I am planning to make a tutorial about Helmholtz equation applied on computational acoustics and a comparison with commercial software (Actran, Comsol, Simcenter 3D, etc.)

5

Hi kamensky,
How can I check which quadrature rule I am using? What is the default value for standard lagrange elements 1st and 2nd order?

6

The quadrature rule depends on the variational form, and is estimated by ufl if not explicitly supplied. See: How many integration points are included in the tetrahedral grid? - #2 by dokken for more details

1 Like
7

Hi dokken.
Are there news about Point Source reimplementation?
Thank you

8

The parallel bit of it required some more thought than I first expected.

However, if you want something that works in serial, I can give you a simply implementation for that.

1 Like
9

Yes! It would be wonderful!

10

Hi Everyone,

Wanted to ask if PointSource has been completed or is nearing completion anytime soon? I suppose I could go back to fenics but it would have been nice to use fenicsx as I have a preference over hexahedral meshes at the moment (imported from gmsh which I understand is not supported by fenics). Greatly appreciate it.

11

Hi everyone,
I am also interested in this topic. Currently, I am in the process of migrating some scripts that utilize point sources from FEniCS to FEniCSx. I tried to add them with methods described in the topic Dirac delta distribution (dolfinx) but so far without success. Reimplementation of PointSource would be very helpful.

12

Hi dokken,

Would it be possible to get a copy of the simple implementation as well? Iā€™d greatly appreciate it.

13

Hi Everyone,

I think that I made a progress. By combining information from Dirac delta distribution (dolfinx) with Cahn-Hilliard equation demo I realized that u function in Dokkenā€™s example is not TrialFunction as it is in many demos, but a separate function.

I attached a demo code below. Letā€™s assume that you have a square 10x10 domain, and want to put two sources of intensities 1 and -1 on coordinates (2.5, 5) and (7.5, 5).

import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (10.0, 10.0)), n=(100, 100),
                            cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))

facets = mesh.locate_entities_boundary(msh, dim=1, marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 10.0)),
                                                                                  np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 10.0))))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)

f = fem.Function(V)
dofs = fem.locate_dofs_geometrical(V,  lambda x: np.isclose(x.T, [2.5, 5, 0]).all(axis=1))
f.x.array[dofs] = 1
dofs = fem.locate_dofs_geometrical(V,  lambda x: np.isclose(x.T, [7.5, 5, 0]).all(axis=1))
f.x.array[dofs] = -1

a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

try:
    import pyvista
    pyvista.start_xvfb()
    cells, types, x = plot.create_vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(cells, types, x)
    grid.point_data["u"] = uh.x.array.real
    grid.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=False)
    warped = grid.warp_by_scalar()
    plotter.add_mesh(warped)
    plotter.show(cpos='xy')
except ModuleNotFoundError:
    print("'pyvista' is required to visualise the solution")
    print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")

Result of modelling:

Result
Result1024Ɨ768 34.6 KB

Please note that this method works only if the location of point sources aligns with the mesh.
Right now I am not entirely sure if everything works correctly but the results look promising. I would be very grateful for any confirmation or corrections.

1 Like
14

Hi eMWu,

Your implementation is in fact how I try it currently and I get similar results. I am using a structured hexagonal grid for a cube with the source point at the origin. Thus, my source point coincidences with a grid point as well. The problem for me is that my numerical solution from dolfinx is off from the analytical solution by a scalar factor. The profile is the same - a diffusive solution with smooth derivatives (away from the source point). The scalar factor difference is orders of magnitude so it is significant. The numerical solution also changes significantly depending on the size of the mesh. So I think the problem is how the Dirac delta function is being interpolated onto the grid and function space, something that PointSource class in Fenics would have taken care of. For example in finite difference, we simply need to divide the source term by the size of the cell.

Have you tried comparing your results to analytical solutions? Or looking at the mesh sensitivity of the solution?

15

Hi, yes after posting I compared numerical and analytical solutions for a point source of current located in a homogenous medium of constant conductivity. For the numerical solution, I used equations for the radially symmetric model.
Unfortunately, you are right. Assuming that I didnā€™t make any other mistakes numerical solutions underestimate the potential around the source by several orders of magnitude. It is also changing with the mesh size. Curiously, when the source is added with the method described above the difference between analytical and numerical solutions is greater when a finer mesh Is present around the source. When the source is added as an approximation of the Dirac delta the mesh size is not affecting the results as much.

Result
Result896Ɨ876 81.2 KB

16

When you add it as a degree of freedom of a first order space, you are approximating the dirac delta function as a piecewise linear function going from one at the dof to zero at all other dofs of the cell. Thus the result will be highly dependent on the mesh resolution.

If you use a dirac delta function, it will be evaluated at quadrature points. These are also mesh dependent, but not as prominent.

17

So I take it that the best way in dolfinx at the moment is to use the approximation? Or would there be an implementation of it that you mentioned earlier in the thread that could be a better solution @dokken?

Also, is there a good rule to follow for the value of ā€˜aā€™ in the approximation that would depend on the mesh size?

18

Ive got alot of things going on at the moment, so i havenā€™t had time to get pointsources to work nicely in parallel. Iā€™ll hopefully have more time in a few weeks to work on it.

19

Thank you again for your help. I was talking more about the implementation you mentioned you already had working in serial? But I understand you are very busy.

20

Hi, thank you for the answers. I am not in rush with conversion to fenicsx so probably I will just wait for the implementation of point sources.