Trying to solve a homogeneous Poisson equation with point sources

General
1

Hello,

I am trying to solve a homogeneous (insulated) Poisson equation similar to the example real_function_space.py on scifem . 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.

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

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:

# 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.

2

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:

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)
1 Like
3

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.