Why is the Laplacian of my computed solution to Laplace's Equation nonzero?

[fea]

I have a harmonic vector \textbf{V} = (x^2-y^2)\hat{\textbf{i}} + (x^2-y^2)\hat{\textbf{j}} in a square shaped domain. I used this harmonic vector to formulate the boundary conditions for a scalar potential by computing the flux at each point
\phi(s) = \phi(s_0) - \int_{s_0}^s \textbf{V}\cdot \hat{\textbf{n}}ds
and I want my scalar potential to be harmonic as well i.e. \nabla^2 \phi = 0. Unfortunately, the Laplacian of this solution fails to satisfy the right-hand side of the Laplace Equation almost everywhere inside the domain. Below is a plot of where Laplace Equation being satisfied (\nabla^2 \psi = 0) was true, and where it was false (\nabla^2 \phi \ne 0)

Here is my code:

from mpi4py import MPI
import dolfinx as df
import ufl
import numpy as np
from dolfinx.fem.petsc import LinearProblem
from petsc4py import PETSc
import pyvista as pv

class bc_setter:
    def __init__(self, V_h):
        self.V_h = V_h
        self.axis = 0  # initial axis of integration
        # initial x_c (will be adjusted in cumul_int_expr for each quadrature point)
        self.x_c = df.fem.Constant(mesh, PETSc.ScalarType(0.0))
        self.n = ufl.FacetNormal(mesh)  

        # define all boundary facets for integration (same markers you used)
        self.boundaries = [
            (1, lambda x: np.isclose(x[1], -Ly)),  # Southern boundary
            (2, lambda x: np.isclose(x[0], Lx)),   # Eastern boundary
            (3, lambda x: np.isclose(x[1], Ly)),   # Northern boundary
            (4, lambda x: np.isclose(x[0], -Lx)),  # Western boundary
        ]
        # build meshtags for boundary markers (same as your previous logic)
        facet_indices, facet_markers = [], []
        for marker, locator in self.boundaries:
            facets = df.mesh.locate_entities(mesh, fdim, locator)
            facet_indices.append(facets)
            facet_markers.append(np.full_like(facets, marker))

        facet_indices = np.hstack(facet_indices).astype(np.int32)
        facet_markers = np.hstack(facet_markers).astype(np.int32)
        sorted_facets = np.argsort(facet_indices)
        self.facet_tag = df.mesh.meshtags(
            mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
        )
        self.ds = ufl.Measure("ds", domain=mesh, subdomain_data=self.facet_tag)

        # initialize a placeholder compiled form;
        # note: this will be overwritten in set_psi_bcs() for each boundary segment
        self.compiled_form = df.fem.form(
            (ufl.dot(self.n, self.V_h))
            * ufl.conditional(self.x_c < 0, 1, 0)
            * ufl.ds
        )

    def cumul_int_expr(self, x):
        """Function used by interpolate: computes cumulative integral up to x[self.axis]."""
        out_arr = np.zeros(x.shape[1], dtype=np.float64)
        for i in range(x.shape[1]):
            self.x_c.value = x[self.axis][i] # set the cutoff location for the conditional in the form
            val = df.fem.assemble_scalar(self.compiled_form)
            val = mesh.comm.allreduce(val, op=MPI.SUM)
            out_arr[i] = val
        return out_arr
    #---- Create boundary conditions ----
    def set_psi_bcs(self, V_hn = None):
        """
        Build Dirichlet BC objects for psi_h on the boundary DOFs
        using cumulative integrals of V_h·n.
        """
        if V_hn is None:
            V_hn = ufl.dot(self.n, self.V_h)
        else:
            pass
        
        psi_h_bcs = []
        bf = []  # Diagnostic: raw boundary fluxes 

        for boundary in self.boundaries:
            boundary_marker = boundary[0]
            # choose axis and param s (x or y) for this boundary
            if boundary_marker % 2 == 0:
                self.axis = 1  # y-axis (east/west boundaries param is y)
                s = ufl.SpatialCoordinate(mesh)[1]
            else:
                self.axis = 0  # x-axis (north/south boundaries param is x)
                s = ufl.SpatialCoordinate(mesh)[0]

            # For the cumulative integral, choose integration direction convention and set compiled_form
            if boundary_marker <= 2:  # positive integration direction
                endpoint_index = -1
                # NOTE: subtract mean_flux_density here to enforce compatibility
                self.compiled_form = df.fem.form(
                    -(V_hn)
                    * ufl.conditional(s < self.x_c, 1, 0)
                    * self.ds(boundary_marker)
                )
            else:  # negative integration direction
                endpoint_index = 0
                self.compiled_form = df.fem.form(
                    -(V_hn)
                    * ufl.conditional(s < self.x_c, 0, 1)
                    * self.ds(boundary_marker)
                )
            # evaluate psi_h(s_0)
            if boundary_marker == 1:
                psi_h_s0 = 0.0
            else:
                # find boundary DOFs for S on this facet group
                boundary_facets = self.facet_tag.find(boundary_marker)
                boundary_dofs = df.fem.locate_dofs_topological(S, fdim, boundary_facets)
                # sort DOFs along the parametric axis to pick endpoint
                coords = S.tabulate_dof_coordinates()[boundary_dofs]
                sorted_indices = np.argsort(coords[:, self.axis])
                psi_h_s0 = psi_h_s.x.array[boundary_dofs][sorted_indices][endpoint_index]

            # Build psi_h_s via interpolate of cumul_int_expr
            psi_h_s = df.fem.Function(S)
            psi_h_s.interpolate(self.cumul_int_expr)
            psi_h_s.x.array[:] = psi_h_s0 + psi_h_s.x.array[:]

            # ---- DIAGNOSTIC: check flux of V_h over each boundary ----
            compiled_form_n = df.fem.form(V_hn * self.ds(boundary_marker))
            raw_flux = mesh.comm.allreduce(df.fem.assemble_scalar(compiled_form_n), op=MPI.SUM)
            bf.append(raw_flux)

            # Decide interpolation DOFs for Dirichlet BC 
            if boundary_marker<3:
                interpolation_facets = self.facet_tag.find(boundary_marker)
                interpolation_dofs = df.fem.locate_dofs_topological(S, fdim, interpolation_facets)    
            else: 
                integration_facets = self.facet_tag.find(boundary_marker)
                integration_dofs = df.fem.locate_dofs_topological(S, fdim, integration_facets)
                sorted_indices = np.argsort(S.tabulate_dof_coordinates()[integration_dofs][:, 1])
                dof_to_remove = integration_dofs[sorted_indices][0]
                A = integration_dofs.tolist(); A.remove(dof_to_remove)
                interpolation_dofs = np.array(A, dtype=np.int32)
            psi_h_bcs.append(df.fem.dirichletbc(psi_h_s, interpolation_dofs))
            
        # Print diagnostic
        print("raw fluxes of V_h along 4 boundaries = ", bf)
        return psi_h_bcs

# ---- Diagnostic functions ----
def compute_flux(F):
    # F can be a Function in V space or a UFL vector expression
    form = df.fem.form(ufl.dot(n, F) * ufl.ds)
    val = df.fem.assemble_scalar(form)
    return mesh.comm.allreduce(val, op=MPI.SUM)
def compute_circulation(F):
    form = df.fem.form(ufl.dot(t, F) * ufl.ds)
    val = df.fem.assemble_scalar(form)
    return mesh.comm.allreduce(val, op=MPI.SUM)

Nx,Ny = 60,60 # Number of points in x and y
Lx, Ly = 1,1
mesh = df.mesh.create_rectangle(MPI.COMM_WORLD, [[-Lx,-Ly],[Lx,Ly]] ,[Nx, Ny])
tdim = mesh.topology.dim
fdim = tdim-1
mesh.topology.create_connectivity(tdim-1, tdim)

# Define normal and tangential unit vectors
n = ufl.FacetNormal(mesh)
t = ufl.perp(n)

S = df.fem.functionspace(mesh, ("CG", 2))
V = df.fem.functionspace(mesh, ("CG", 2, (2,))) #order CG satisfies all constraints

# Exact harmonic vector
Vex = df.fem.Function(V)
k1, k2 = 1, 1
Vex.interpolate(lambda x: (k1*(x[0]**2-x[1]**2),k2*(x[0]**2 - x[1]**2)))

#---- Diagnostic: confirm vector is harmonic ----
print("flux of harmonic vector (should be ~0) = ",compute_flux(Vex))
print("circulation of harmonic vector (should be ~0) = ",compute_circulation(Vex))

# Extract boundary conditions from harmonic vector
bc_obj = bc_setter(Vex)
bc_stream = bc_obj.set_psi_bcs()

# ---- Laplace Equation ----
u = ufl.TrialFunction(S)
v = ufl.TestFunction(S) 
a = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx # a(u,v)
f = df.fem.Constant(mesh,df.default_scalar_type(0))
L_stream = -f*v*ufl.dx 
problem = LinearProblem(a=a, L=L_stream, bcs = bc_stream, petsc_options_prefix="Laplace",petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) 
psi_h = problem.solve()

# -------------------------------------
# ---- Diagnostics: Test solution -----
# -------------------------------------

#---- Diagnostic: Test Green's Theorem ----
print("\nCheck Green's Theorem\n---------------------")
print("int int nabla^2 u = ",df.fem.assemble_scalar(df.fem.form(ufl.div(ufl.grad(psi_h)) * ufl.dx)))
print("oint du/dn ds = ",compute_flux(ufl.grad(psi_h)))
diff = df.fem.assemble_scalar(df.fem.form(ufl.div(ufl.grad(psi_h)) * ufl.dx)) - compute_flux(ufl.grad(psi_h))
print("difference = ",diff)
if np.isclose(diff, 0):
    print("---------------------\nGreen's Theorem passed\n")
else:
    print("Green's Theorem failed\n")

# ---- Check: satisfies Laplace Equation ----
lap_psi = df.fem.Function(S)
lap_expr = df.fem.Expression(ufl.div(ufl.grad(psi_h)), S.element.interpolation_points)
lap_psi.interpolate(lap_expr)
bool_array = np.isclose(lap_psi.x.array[:], 0, atol= 10E-6)

print("Check Laplace's Equation\n---------------------")
lap_psi = df.fem.Function(S)
lap_expr = df.fem.Expression(ufl.div(ufl.grad(psi_h)), S.element.interpolation_points)
lap_psi.interpolate(lap_expr)
if np.mean(np.all(np.isclose(lap_psi.x.array[:], 0))) > .95: # If roughly 95% of points satisfy Laplace's Equation, it passes
    print("Laplace's Equation passed")
else:
    print(f"Laplace's Equation failed. {100*(1-np.mean(np.isclose(lap_psi.x.array[:], 0, atol = 10E-8))) :.2f}% of points failed to satisfy Laplace's Equation.\nShowing plot where Laplace's Equation passed and failed:")
    grid = pv.UnstructuredGrid(*df.plot.vtk_mesh(S))
    grid["u"] = bool_array
    plotter = pv.Plotter()
    plotter.add_mesh(grid, cmap = "bwr")
    plotter.view_xy()
    plotter.show_grid()
    plotter.show()

Are my Dirichlet boundary conditions at fault; is there some sort of condition they must first satisfy for the solver to satisfy the PDE?

I would also like to know if changing my solver configuration in any way might help.

[Stein]

Finite element schemes to not aim to satisfy the PDE pointwise. Their weak-form origin means the PDE is only satisfied in a “weighted average sense”. So, you would not actually expect your solution to satisfy \Delta \phi^h = 0 (if it would, you would have found the exact solution to the PDE, not the discrete one). Of course, with mesh refinement, it would get closer and closer.

[fea]

Oh okay. Raising the order did in fact cause it to converge to approximately satisfying the Laplace’s Equation (in a majority of the domain)
3rd-order:

4th-order:

5th-order:

5th-order is pretty time-expensive though. I might try to see if 4th-order is good enough, and even then it would still take a while for postprocessing. For faster postprocessing with a more accurate solution, would interpolating this solution onto a lower order space like for example, moving my 5th-order solution to a 2nd-order space, work? Or would it lose its accuracy doing that?

[Stein]

I can only imagine that you would indeed see a significant loss of accuracy. There is nothing in the interpolation logic that emphasizes the importance of matching \Delta \phi^h.

But I think the more important question to ask yourself is why you’re so interested in matching \Delta \phi^h = 0 pointwise, and if you can avoid that necessity. “Collocation” methods do satisfy that requirement at certain points (the collocation points), but no discrete method will give you that at every point. Unless… you’ve retrieved the exact solution.