Poisson equation with non-zero right hand side and non-zero neumann boundary condition in one section of an edge and zero neumann boundary condition everywhere else

[abhishek-appana]

Hello,

I was using FeNICSx earlier and decided to switch to scikit-fem recently. I am trying to implement the following problem based on the examples I have seen till now. But I am not sure where I am going wrong. Would be glad if someone could help.

The problem is :
$-\Delta u = f$
$f$ is independent of $u$

The geometry and boundaries are defined as:

The boundaries 1, 2, 3, 5, and 6 have zero neumann boundary conditions, that is, $\partial u/\partial n = 0$, where $n$ is the direction normal to the edge
Boundary 4 has a non-zero neumann boundary condition, that is, $\partial u/\partial n = -r*u$, $r$ is independent of $u$

The weak variational formulation that comes from this is:
LHS: $A(u,v) = \int _\Omega (\nabla u . \nabla v) dx + \int _{4} r * u * v ds$, where $\Omega$ is the whole domain and $4$ is the edge which is labelled $4$ in the geometry

RHS: $b(v) = \int _\Omega f*v dx$

The code I wrote for this implementation is as below:

from skfem import *
from skfem.helpers import dot, grad
from skfem.models.poisson import laplace, mass
from math import ceil
import numpy as np

length = 0.193
height = 0.195
ncells = 50

I = 2
Nlayers = 20
Acs = 0.193*0.195
Sigma = 3.5e7
delta = 20e-6
Acst = 0.018*delta

f = I/(Nlayers*Acs*Sigma*delta)
r = -1*I/(Nlayers*Acst*Sigma)

mesh = (MeshTri
        .init_tensor(np.linspace(0, length, ncells),
                     np.linspace(0, height, ncells))
        .with_boundaries({'1': lambda x: np.isclose(x[1],0),
                          '2': lambda x: np.isclose(x[0],0),
                          '3': lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0) & np.less_equal(x[0],0.119),
                          '4': lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0.119) & np.less_equal(x[0],0.169),
                          '5': lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0.169) & np.less_equal(x[0],0.193),
                          '6': lambda x: np.isclose(x[0],0.193)}))

e = ElementTriP2()
basis = {
    'cell': Basis(mesh, e),
    'tab': FacetBasis(mesh, e, facets=mesh.boundaries['4'])}

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))

@BilinearForm
def tab(u,v,_):
    return r*u*v

A = (asm(laplace, basis['cell']) + asm(tab, basis['tab']))

@LinearForm
def rhs(v, _):
    return f*v

b = asm(rhs, basis['cell'])

D = basis['cell'].get_dofs()
I = basis['cell'].complement_dofs(D)

u = solve(*condense(A, b, I=I))

def visualize():
    from skfem.visuals.matplotlib import plot, show
    return plot(basis['cell'], u, shading='gouraud', colorbar=True)

if __name__ == '__main__':
    visualize().show()

When I run this, I get the output:

I am expecting the output:

Thank you for the help!

[kinnala]

I think you don't need condense() now because it is used for essential boundary conditions. So just solve should be enough.

[kinnala]

If you can post the corresponding FEniCS code then I can check if there is any obvious difference?

[kinnala]

Alright, now I see this could also be related to the same issue: how the facets are defined in scikit-fem vs in FEniCS when indicator function is used.

Since we are close to the blue point, where y-gradient changes from one to another due to the change in the BC, I can see that the details will determine the actual value of the gradient. To the left from the blue point the gradient is zero. To the right there is a Robin condition and nonzero gradient. I will make some tests.

[kinnala]

The closest match (to FEniCS) I am able to get with this small change:

                              '4': lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0.119 + length/ncells) & np.less_equal(x[0],0.169),

It produces the following $u$:

The gradient values are:
dx: [-0.00296201], dy: [-0.00023504]

I observed that the y-gradient is highly dependent on the details on how this integration over "4" is performed. If I wanted the best possible results (by "best" I do not mean closest to FEniCS but closest to the true solution of the underlying PDE), I would probably design a mesh which has these endpoints of "4" as vertices, to make sure that the integration over "4" is performed exactly. However, it is up to you whether you want more accurate result or not. All approaches discussed so far do not seem to lead to an accurate result corresponding to the underlying PDE.

[abhishek-appana]

Hello @kinnala,

Apologies for the delayed response!

Firstly, here is the FEniCSx code.

# %%
from dolfinx import default_scalar_type, geometry
import dolfinx
from dolfinx.fem import (Constant, Function, FunctionSpace, form, locate_dofs_topological, Expression)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, locate_entities, meshtags, CellType
from dolfinx.plot import vtk_mesh

from mpi4py import MPI
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction, dx, grad, inner, lhs, rhs, as_vector)

from petsc4py import PETSc
from mpi4py import MPI

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pyvista

# %% [markdown]
# Creating a rectangular geometry with a triangular mesh

# %%
nx, ny = 50, 50
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]),np.array([0.193,0.195])], [nx, ny], CellType.triangle)
V = FunctionSpace(mesh, ("Lagrange",1))

# %% [markdown]
# Defining the spatial coordinates

# %%
x = SpatialCoordinate(mesh)

# %%
I = Constant(mesh,default_scalar_type(2))
Nlayers = Constant(mesh,default_scalar_type(20))
Acs = Constant(mesh,default_scalar_type(0.193*0.195))
Sigma = Constant(mesh,default_scalar_type(3.5e7))
delta = Constant(mesh,default_scalar_type(20e-6))
Acst = Constant(mesh,default_scalar_type(0.018*delta))
f = Constant(mesh,default_scalar_type(I/(Nlayers*Acs*Sigma*delta)))
r = Constant(mesh,default_scalar_type(-1*I/(Nlayers*Acst*Sigma)))

# %% [markdown]
# Creating boundary definitions

# %%
boundaries = [(1, lambda x: np.isclose(x[1],0)),
              (2, lambda x: np.isclose(x[0],0)),
              (3, lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0) & np.less_equal(x[0],0.119)),
              (4, lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0.119) & np.less_equal(x[0],0.169)),
              (5, lambda x: np.isclose(x[1],0.195) & np.greater_equal(x[0],0.169) & np.less_equal(x[0],0.193)),
              (6, lambda x: np.isclose(x[0],0.193))]

# %% [markdown]
# Creating meshtags

# %%
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = 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)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

# %% [markdown]
# Debugging boundary conditions

# %%
with XDMFFile(mesh.comm, "facet_tags_CB.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tag, mesh.geometry)

# %%
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# %% [markdown]
# Initializing variational formulation

# %%
u, v = TrialFunction(V), TestFunction(V)
# a = form(inner(grad(u), grad(v)) * dx + r * u * v * ds(4))
# L = form(f * v * dx)
F = (inner(grad(u), grad(v)) * dx + r * u * v * ds(4) - f * v * dx)
a = lhs(F)
L = rhs(F)

# %% [markdown]
# Defining solver and solving the problem

# %%
CB_solve = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = CB_solve.solve()

# %% [markdown]
# Visualization

# %%
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")

plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("Charge balance.png")

# %% [markdown]
# Calculating gradients at various points

# %%
x_grad = grad(uh)[0]
y_grad = grad(uh)[1]
W = FunctionSpace(mesh, ("DQ", 1))
xgrad_expr = Expression(x_grad, W.element.interpolation_points())
ygrad_expr = Expression(y_grad, W.element.interpolation_points())
xgrad = Function(W)
xgrad.interpolate(xgrad_expr)
ygrad = Function(W)
ygrad.interpolate(ygrad_expr)

# %%
points = np.array([[0.119],  [0.1945], [0]])

# %%
bb_tree = geometry.bb_tree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])

# %%
points_on_proc = np.array(points_on_proc, dtype=np.float64)
xgrad_values = xgrad.eval(points_on_proc, cells)
ygrad_values = ygrad.eval(points_on_proc, cells)
print("Gradient in X = ", xgrad_values, "; Gradient in Y = ", ygrad_values)

[abhishek-appana]

Thank you so much for trying various meshes and providing possible solutions @kinnala.

I need to increase the accuracy and get the closest to the true solution of PDE. But I cannot change the position of edge 4, like you have suggested.

How can I create a mesh to ensure that the endpoints of 4 match the vertices of an element? And how to create a mesh with more elements close to this area (let's say, in the right top half of the geometry) and fewer elements away from it?