Laplace on 3d sphere/ellipsoid - exterior problem

variational formulation
1

I’m working on solving the Laplace equation on a sphere with simple Dirichlet (constant) and Neumann(0) boundary conditions. Analytically the second boundary condition is a constant potential at infinity, I didnt know how to enforce this condition so I tried to use the 1/r decay of the potential at large distances from the sphere.

Currently, I use two coincide spheres: an inner sphere with the correct boundary conditions and a larger sphere (with a 1:10 radii ratio) to simulate the 1/r decay. While this approach works for various tests I performed, I’m unsure if it’s correct and efficient solution, also different geometries like an ellipsoid might need a different approach.

Has anyone dealt with a similar setup or knows of a more optimized method to handle the 1/r decay (/ constant at infinity)?
I’d appreciate any insights!

import dolfinx.io.gmshio
import numpy as np
from mpi4py import MPI
import gmsh
import ufl
from dolfinx import fem, io, plot, mesh
from dolfinx.fem import dirichletbc, locate_dofs_geometrical, FunctionSpace, Constant
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunction, TestFunction, grad, inner, ds, dx
import pyvista
from dolfinx.plot import vtk_mesh as create_vtk_mesh
from utils import *

shell_volume_tag = 1
inner_boundary_tag = 2
outer_boundary_tag = 3
REFINEMENT_RADIUS = 0.1

def generate_spherical_shell_mesh(a, R, inner_mesh_size, outer_mesh_size):

    return msh, cell_tags, facet_tags


# Mesh parameters
a = 1.0  # Cell radius
R = 10.0  # Outer boundary radius
inner_mesh_size = 0.05  # Refined mesh near receptor
outer_mesh_size = 1.0
phi_0 = 1.0

gmsh.initialize()
gmsh.model.add("spherical_shell")
# Create the inner sphere (volume)
inner_sphere = gmsh.model.occ.addSphere(0, 0, 0, a)
gmsh.model.occ.synchronize()
# Create the outer sphere (volume)
outer_sphere = gmsh.model.occ.addSphere(0, 0, 0, R)
gmsh.model.occ.synchronize()
# Create the spherical shell by cutting the inner sphere from the outer sphere
shell = gmsh.model.occ.cut([(3, outer_sphere)], [(3, inner_sphere)])
gmsh.model.occ.synchronize()
# Get the volume tag of the spherical shell
shell_volume = shell[0][0][1]
# Get the surfaces of the spherical shell
shell_surfaces = gmsh.model.occ.getEntities(dim=2)
gmsh.model.occ.synchronize()
# Add physical groups
gmsh.model.addPhysicalGroup(3, [shell_volume], tag=shell_volume_tag)  # Spherical shell volume
gmsh.model.setPhysicalName(3, 1, "ShellVolume")
gmsh.model.addPhysicalGroup(2, [inner_sphere], tag=inner_boundary_tag)  # Inner boundary
gmsh.model.setPhysicalName(2, 1, "InnerBoundary")
gmsh.model.addPhysicalGroup(2, [outer_sphere], tag=outer_boundary_tag)  # Outer boundary
gmsh.model.setPhysicalName(2, 2, "OuterBoundary")

gmsh.model.occ.synchronize()
# Set mesh sizes on the surfaces
gmsh.model.mesh.setSize(gmsh.model.getBoundary([(2, inner_sphere)], recursive=True),inner_mesh_size)
gmsh.model.mesh.setSize(gmsh.model.getBoundary([(2, outer_sphere)], recursive=True),outer_mesh_size)

# Mesh generation
gmsh.model.mesh.generate(3)
# Convert to DOLFINx mesh
msh, cell_tags, facet_tags = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()



V = fem.functionspace(msh, ('CG', 2))
# Define measures with facet_tags
inner_tag = inner_boundary_tag
outer_tag = outer_boundary_tag
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags)
dx = ufl.dx(domain=msh)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Locate dofs on the tagged surface using geometrical condition
inner_dofs = fem.locate_dofs_topological(V, 2, facet_tags.find(inner_boundary_tag))
outer_dofs = fem.locate_dofs_topological(V, 2, facet_tags.find(outer_boundary_tag))

# Define Dirichlet BC on inner sphere
bc_inner = fem.dirichletbc(value=phi_0, dofs=inner_dofs, V=V)

# Define Neumann condition (zero flux) on the inner boundary
# Neumann conditions with zero flux should be naturally handled, but we include it explicitly for clarity
g_N = Constant(msh, 0.0)

# Define the variational forms
a_form = inner(grad(u), grad(v)) * dx + (1 / R) * u * v * ds(outer_tag)
L_form = g_N * v * ds(inner_tag) + Constant(msh, 0.0) * v * dx

# Set up and solve the linear problem
problem = LinearProblem(
    a_form,
    L_form,
    bcs=[bc_inner],
    petsc_options={
        "ksp_type": "cg",
        "pc_type": "hypre",
        "ksp_rtol": 1e-9,
        "ksp_max_it": 1000,
        "ksp_monitor": None
    }
)
uh = problem.solve()

# flux calculation
n = ufl.FacetNormal(msh)
flux_form_inner = inner(grad(uh), n) * ds(inner_tag)
total_flux_inside = fem.assemble_scalar(fem.form(flux_form_inner))
MPI.COMM_WORLD.barrier()
# Compute the analytical flux for comparison
flux_inner_exact = 4 * np.pi * phi_0 * a
print(f"Numerical flux over the inner sphere: {total_flux_inside:.6f}")
print(f"Analytical flux over the inner sphere: {flux_inner_exact:.6f}")

the results looks good:

Numerical flux over the inner sphere: 12.507325
Analytical flux over the inner sphere: 12.566371
2

Hello,

This sounds like the kind of problem that the boundary element method (BEM) might be better for, as BEM can naturally handle infinite domains.

Taking a look at the Bempp tutorials might be a good place to start.

2 Likes