# Trying to solve a homogeneous Poisson equation with point sources

**URL:** <https://fenicsproject.discourse.group/t/trying-to-solve-a-homogeneous-poisson-equation-with-point-sources/16896>\
**Category:** General\
**Tags:** dolfinx\
**Created:** [February 1, 2025, 6:29am UTC](https://fenicsproject.discourse.group/t/trying-to-solve-a-homogeneous-poisson-equation-with-point-sources/16896 "2025-02-01T06:29:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![prjohnston](https://avatars.discourse-cdn.com/v4/letter/p/5e9695/32.png) [@prjohnston](https://fenicsproject.discourse.group/u/prjohnston)\
**Post date:** [February 1, 2025, 6:29am UTC](https://fenicsproject.discourse.group/t/trying-to-solve-a-homogeneous-poisson-equation-with-point-sources/16896/1 "2025-02-01T06:29:08Z")

</div>

Hello,

I am trying to solve a homogeneous (insulated) Poisson equation similar to the example [real\_function\_space.py](https://scientificcomputing.github.io/scifem/examples/real_function_space.html) on [scifem](https://scientificcomputing.github.io/scifem/README.html) . In my example the both the `g` and `h` functions are identically zero and the function `f` is a sum of two Dirac delta functions.

I am trying to model my solution approach on the scifem example [point\_source.py](https://scientificcomputing.github.io/scifem/examples/point_source.html).

However, I cannot merge the two approaches together. I must be missing some subtle point.

From the `real_function_space.py` example I have the variational problem

```auto
a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx
L0 = ufl.inner(f, du) * ufl.dx
L1 = ufl.inner(zero, dl) * ufl.dx
a = dolfinx.fem.form([[a00, a01], [a10, None]])
L = dolfinx.fem.form([L0, L1])

```

Next, following the `point_source.py` example, I created the source and sink points as follows:

````auto
# We are now ready to apply the point source.
# First, we create the right hand side vector for the problem, and initialize it as zero.

bp = dolfinx.fem.Function(V)
bp.x.array[:] = 0

# Secondly we define the point sources we want to apply.
#
# ```{warning}
# Note that if running in parallel, we only supply the points on a single process.
# The other process gets an empty array.
# ```

geom_dtype = mesh.geometry.x.dtype

# firstly the source point

if mesh.comm.rank == 0:
    points = np.array([[0.475, 0.5, 0]], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)

# Next, we create the point source object and apply it to the right hand side vector.

gamma = 1.
point_source = scifem.PointSource(V, points, magnitude=gamma)
point_source.apply_to_vector(bp)

# now the sink point

if mesh.comm.rank == 0:
    points = np.array([[0.525, 0.5, 0]], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)

gamma = -1.
point_source = scifem.PointSource(V, points, magnitude=gamma)
point_source.apply_to_vector(bp)

````

My problem is that I can’t work out how to incorporate the information in `bp` into the form `L`.

I did manage to work out that I could use this definition of `bp` as `f` in the expression for `L0` above and everything works fine. However, I didn’t think that this was the way the `point_source` function was supposed to work.

I would be grateful if someone could provide me with some insights as to the points I am missing or correct the mistakes that I have made.

BTW, I am using fenicsx 0.9.0 and scifem 0.3.0.

Thank you very much for taking the time read this long message. Any insights would be much appreciated.

Peter.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 1, 2025, 3:51pm UTC](https://fenicsproject.discourse.group/t/trying-to-solve-a-homogeneous-poisson-equation-with-point-sources/16896/2 "2025-02-01T15:51:25Z")

</div>

You can spot this in the original real space demo, where one does `scatter_local_vectors`. By using this command, you can insert `bp` into the assembled RHS.  
Here is a minimal reproducible (working) example:

```python
from packaging.version import Version
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors
import dolfinx.fem.petsc

import numpy as np
from scifem import create_real_functionspace, assemble_scalar
import ufl

M = 40
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

def u_exact(x):
    return 0.3 * x[1] ** 2 + ufl.sin(2 * ufl.pi * x[0])

x = ufl.SpatialCoordinate(mesh)
n = ufl.FacetNormal(mesh)
g = dolfinx.fem.Constant(mesh, 0.0)
f = dolfinx.fem.Constant(mesh, 0.0)
h = assemble_scalar(u_exact(x) * ufl.dx)

# We then create the Lagrange multiplier space

R = create_real_functionspace(mesh)

# Next, we can create a mixed-function space for our problem

if dolfinx. __version__ == "0.8.0":
    u = ufl.TrialFunction(V)
    lmbda = ufl.TrialFunction(R)
    du = ufl.TestFunction(V)
    dl = ufl.TestFunction(R)
elif Version(dolfinx. __version__ ) >= Version("0.9.0.0"):
    W = ufl.MixedFunctionSpace(V, R)
    u, lmbda = ufl.TrialFunctions(W)
    du, dl = ufl.TestFunctions(W)
else:
    raise RuntimeError("Unsupported version of dolfinx")

# We can now define the variational problem

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))

a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx
L0 = ufl.inner(f, du) * ufl.dx + ufl.inner(g, du) * ufl.ds
L1 = ufl.inner(zero, dl) * ufl.dx

a = dolfinx.fem.form([[a00, a01], [a10, None]])
L = dolfinx.fem.form([L0, L1])

# Note that we have defined the variational form in a block form, and
# that we have not included $h$ in the variational form. We will enforce this
# once we have assembled the right hand side vector.

# We can now assemble the matrix and vector

A = dolfinx.fem.petsc.assemble_matrix_block(a)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector_block(L, a, bcs=[])

# Next, we modify the second part of the block to contain `h`
# We start by enforcing the multiplier constraint $h$ by modifying the right hand side vector

if dolfinx. __version__ == "0.8.0":
    maps = [
        (V.dofmap.index_map, V.dofmap.index_map_bs),
        (R.dofmap.index_map, R.dofmap.index_map_bs),
    ]
elif Version(dolfinx. __version__ ) >= Version("0.9.0.0"):
    maps = [(Wi.dofmap.index_map, Wi.dofmap.index_map_bs) for Wi in W.ufl_sub_spaces()]

b_local = get_local_vectors(b, maps)
b_local[1][:] = 0

# Next, we create the point source object and apply it to the right hand side vector.

import scifem

b_poisson = dolfinx.fem.Function(V)

geom_dtype = mesh.geometry.x.dtype
if mesh.comm.rank == 0:
    points = np.array([[0.68, 0.362, 0]], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)
point_source = scifem.PointSource(V, points, magnitude=1)
point_source.apply_to_vector(b_poisson)

if mesh.comm.rank == 0:
    points = np.array([[0.14, 0.213, 0]], dtype=geom_dtype)
else:
    points = np.zeros((0, 3), dtype=geom_dtype)
point_source = scifem.PointSource(V, points, magnitude=-1)
point_source.apply_to_vector(b_poisson)

b_local[0][:] = b_poisson.x.array
scatter_local_vectors(
    b,
    b_local,
    maps,
)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# We can now solve the linear system

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")

xh = dolfinx.fem.petsc.create_vector_block(L)
ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Finally, we extract the solution u from the blocked system and compute the error

uh = dolfinx.fem.Function(V, name="u")
x_local = get_local_vectors(xh, maps)
uh.x.array[: len(x_local[0])] = x_local[0]
uh.x.scatter_forward()

with dolfinx.io.VTXWriter(mesh.comm, "uh.bp", [uh]) as bp:
    bp.write(0.0)

```

---

<div class="post-metadata">

**Author:** ![prjohnston](https://avatars.discourse-cdn.com/v4/letter/p/5e9695/32.png) [@prjohnston](https://fenicsproject.discourse.group/u/prjohnston)\
**Post date:** [February 2, 2025, 7:53am UTC](https://fenicsproject.discourse.group/t/trying-to-solve-a-homogeneous-poisson-equation-with-point-sources/16896/3 "2025-02-02T07:53:18Z")

</div>

Thank you so much for that. I am not sure that I would have worked this out for myself but I think I understand thing better now.