Robin BC does not agree with exact solution for 1D diffusion

General
1

Hi there,

I try to model simple 1D diffusion problem with Robin BC using Dokken’s tutorial using DOLFINX v0.9.0 with conda. The Dirichlet and Neumann BCs are all OK but Robin BC does not agree with the analytical solution within the domain although it satisfies the values at the boundaries.

Here is the MWE:

from dolfinx.fem import functionspace
from dolfinx.fem import (Constant,  Function, functionspace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological)
from ufl import TrialFunction, TestFunction, inner, ds, lhs, rhs, dx, Measure, grad
import numpy as np
from dolfinx.fem.petsc import LinearProblem
from dolfinx import default_scalar_type
from dolfinx.mesh import meshtags, locate_entities, create_interval
from mpi4py import MPI


# Geometrical parameters
N = 10
L_total = 1.0

mesh = create_interval(MPI.COMM_WORLD, N, [0.0, L_total])
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
                (2, lambda x: np.isclose(x[0], L_total))]

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(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Define the boundary conditions
u_left = 1.0
D_ij = 1.0 # diffusion coefficient
r_k = -1 # mass transfer coefficient
r = Constant(mesh, default_scalar_type(-D_ij*r_k))
s = 0.3 # concentration outside the surface

degree = 1

V = functionspace(mesh, ("Lagrange", degree))
u, v = TrialFunction(V), TestFunction(V)

f = Constant(mesh, default_scalar_type(0))

F = D_ij * inner(grad(u), grad(v)) * dx - inner(f, v) * dx

ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.x.array[:]=values
            facets = facet_tags.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type
    
boundary_conditions = [
    BoundaryCondition("Dirichlet", 1, u_left),
    BoundaryCondition("Robin", 2, (r, s)),
]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
c_numeric = problem.solve()
x_numeric = np.linspace(0, L_total, len(c_numeric.x.array))


# Analytical solution 
x_analytical = np.linspace(0, L_total, N+1)
c_analytical = (-r_k*s*x_analytical + u_left)/(1-r_k*x_analytical)

import matplotlib.pyplot as plt
plt.plot(x_numeric, c_numeric.x.array, "s-r", label="FEM")
plt.plot(x_analytical, c_analytical, label="trivial")
plt.xlabel("x")
plt.ylabel("concentration (%)")
plt.legend()
plt.show()

giving:

appendix_one_dim_dir_rob
appendix_one_dim_dir_rob1920×1440 108 KB

I would be grateful if you can enlighten me.

2

Are you sure about your analytical solution? With a source f=0, y = (a-b\, c\, x)/(1-c\,x) doesn’t seem to satisfy the ODE internally (maybe I’m overlooking something).

3

I agree with Stein. The source term stemming from that equation is not zero.
It should be: (2 c^2 (a - b))/(c x - 1)^3 (ref: -d^2 (a−bcx)/(1−cx) / dx^2 - Wolfram|Alpha).

4

Thank you for your reply. Let me show my manufactured solution. We solve

\nabla^2 \phi = 0.

At the Robin boundary (x=L), we have

\frac{\partial \phi}{\partial x} (L) = r (\phi - s).

At the Dirichlet boundary, we have

\phi(0) = \phi_{left}.

Let’s assume the solution as

\phi = Ax +B,

If we take the first derivative and then apply Robin BC:

\frac{\partial \phi}{\partial x} (L) = r (\phi - s) = A.

If we apply Dirichlet BC:

\phi(0) = B = \phi_{left}.

Hence, the \phi is:

\phi = r (\phi - s)x + \phi_{left},

after rearranging, we get:

\phi = \frac{-r s x + \phi_{left}}{1-rx}.

Also, changing the analytical solution to the @dokken’s gives

c_analytical = (2*r_k**2*(u_left-s))/(r_k*x_analytical-1)**3

and plotting gives:

MWE
MWE1920×1440 71 KB

Could you please tell me what I am doing wrong here?

5

You start of with postulate

But your final result does not satisfy this structure, so something is off.

should say

\frac{\partial \phi}{\partial x} (L) = r (\phi(L) - s) = r (AL+B - s) = A.

Then you can continue to rearrange.

1 Like
6

You observe here that you have assumed \phi of the form Ax + b, but your final solution is not of that form.

How one chose a manufactured solution goes as follows:
Given: \phi_e = \frac{-r s x + \phi_{left}}{1-rx}
insert this into your PDE to compute what force you expect f to be.
Secondly, insert \phi_e in your Robin condition. to find the appropriate g.

I will illustrate this with a simple example.
Assume u_e(x) = \sin(2\pi *x)
Then f=-\nabla \cdot (D_{ij} \nabla u_e) = 4\pi^2 D_{ij} \sin(2\pi x)
Next you have the relation -D_{ij} \frac{\partial u_e}{\partial n} = \alpha ( u_e - g) .
In 1D, n=(1,) at the right end, and thus
g = u_e +\frac{D_{ij}}{\alpha}\frac{\partial u_e}{\partial x}= \sin(2\pi x) +2\pi \frac{D_ij}{\alpha}\cos(2\pi x)
A simple implementation of this is shown below:

from mpi4py import MPI

from dolfinx.fem import (
    functionspace,
)
import dolfinx

import ufl
from ufl import (
    TrialFunction,
    TestFunction,
    inner,
    ds,
    lhs,
    rhs,
    dx,
    grad,
    SpatialCoordinate,
)
import numpy as np
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_interval


# Geometrical parameters
N = 10
L_total = 1.0

mesh = create_interval(MPI.COMM_WORLD, N, [0.0, L_total])


degree = 1
V = functionspace(mesh, ("Lagrange", degree))
u, v = TrialFunction(V), TestFunction(V)


D_ij = 0.5
x = SpatialCoordinate(mesh)
f_manufactured = 4 * ufl.pi**2 * D_ij * ufl.sin(2 * ufl.pi * x[0])
alpha = 0.1
g_manufactured = ufl.sin(2 * ufl.pi * x[0]) + 2 * ufl.pi * D_ij / alpha * ufl.cos(
    2 * ufl.pi * x[0]
)
F = D_ij * inner(grad(u), grad(v)) * dx - inner(f_manufactured, v) * dx
F += alpha * (u - g_manufactured) * v * ds

# Solve linear variational problem
a = lhs(F)
L = rhs(F)


dofs_right = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))
u_bc = dolfinx.fem.Constant(mesh, 0.0)
bc_right = dolfinx.fem.dirichletbc(u_bc, dofs_right, V)
bcs = [bc_right]
problem = LinearProblem(
    a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
c_numeric = problem.solve()

u_exact = ufl.sin(2 * ufl.pi * x[0])
error_L2 = dolfinx.fem.assemble_scalar(
    dolfinx.fem.form(ufl.inner(c_numeric - u_exact, c_numeric - u_exact) * dx)
)
error_global = mesh.comm.allreduce(error_L2, op=MPI.SUM)
print(f"{N=} L2 error: {np.sqrt(error_global):.3e}")

# Analytical solution
x_numeric = V.tabulate_dof_coordinates()[:, 0]
c_analytical = np.sin(2 * np.pi * x_numeric)

import matplotlib.pyplot as plt

plt.plot(x_numeric, c_numeric.x.array, "s-r", label="FEM")
plt.plot(x_numeric, c_analytical, label="trivial")
plt.xlabel("x")
plt.ylabel("concentration (%)")
plt.legend()
plt.savefig("comparison.png")

which yields

N=10 L2 error: 2.526e-02
N=20 L2 error: 6.357e-03
N=40 L2 error: 1.592e-03
2 Likes
7

Thank you for your replies. I see your point @Stein and thank you for providing that insight. My analytical solution was wrong. It has to be derived like below:

\nabla^2 \phi = 0.

At the Robin boundary (x=L), we have

\frac{\partial \phi}{\partial x} (L) = r (\phi - s).

At the Dirichlet boundary, we have

\phi(0) = \phi_{left}.

Let’s assume the solution as

\phi = Ax +B,

If we apply Dirichlet BC:

\phi(0) = B = \phi_{left}.

If we take the first derivative and then apply Robin BC:

\frac{\partial \phi}{\partial x} (L) = r (\phi(L) - s) = r (AL +\phi_{left} - s) = A.

Then A becomes

A = \frac{r_k(\phi_{left}-s)}{1-r_kL}

hence \phi is

\phi = A = \frac{r_k(\phi_{left}-s)}{1-r_kL}x + \phi_{left}

Putting this into the code

Hence, the \phi is:

\phi = r (\phi - s)x + \phi_{left},

after rearranging, we get:

\phi = \frac{-r s x + \phi_{left}}{1-rx}.

Thank you @dokken for providing a nice example. Please find below the full working code that calculates the L^2 norm:

from dolfinx.fem import functionspace
from dolfinx.fem import (Constant,  Function, functionspace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological)
from ufl import TrialFunction, TestFunction, inner, ds, lhs, rhs, dx, Measure, grad
import numpy as np
from dolfinx.fem.petsc import LinearProblem
from dolfinx import default_scalar_type
from dolfinx.mesh import meshtags, locate_entities, create_interval
from mpi4py import MPI


# Geometrical parameters
N = 10
L_total = 1.0

mesh = create_interval(MPI.COMM_WORLD, N, [0.0, L_total])
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
                (2, lambda x: np.isclose(x[0], L_total))]

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(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Define the boundary conditions
u_left = 1.0
D_ij = 1.0 # diffusion coefficient
r_k = -1 # mass transfer coefficient
r = Constant(mesh, default_scalar_type(-D_ij*r_k))
s = 0.3 # concentration outside the surface

degree = 1

V = functionspace(mesh, ("Lagrange", degree))
u, v = TrialFunction(V), TestFunction(V)

F = D_ij * inner(grad(u), grad(v)) * dx

ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.x.array[:]=values
            facets = facet_tags.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type
    
boundary_conditions = [
    BoundaryCondition("Dirichlet", 1, u_left),
    BoundaryCondition("Robin", 2, (r, s)),
]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
c_numeric = problem.solve()
x_numeric = V.tabulate_dof_coordinates()[:, 0]

x_analytical = np.linspace(0, L_total, N+1)
c_f = lambda x:(r_k*(u_left-s))/(1-r_k*L_total)*x + u_left
u_exact = Function(V)
def c_ufl(x):
    return (r_k*(u_left-s))/(1-r_k*L_total)*x[0] + u_left
u_exact.interpolate(c_ufl)

import matplotlib.pyplot as plt

plt.plot(x_numeric, c_numeric.x.array, "s-r", label="FEM")
plt.plot(x_analytical, c_f(x_analytical), label="trivial")
plt.xlabel("x")
plt.ylabel("concentration (%)")
plt.legend()
plt.savefig("MWE.png", dpi=300)
plt.show()

error_L2 = assemble_scalar(
    form(inner(c_numeric - u_exact, c_numeric - u_exact) * dx)
)
error_global = mesh.comm.allreduce(error_L2, op=MPI.SUM)
print(f"{N=} L2 error: {np.sqrt(error_global):.3e}")

giving
N=10 L2 error: 9.115e-16
and

MWE
MWE1920×1440 93.5 KB

1 Like