Issue with Boundary Conditions in 1D Diffusion PDE Example

[GijsBerkenbosch]

Hi all!

I am trying to solve a 1D diffusion PDE using the py-pde package, but I'm encountering results that differ significantly from my expectations. I suspect the boundary conditions might not be working correctly, or there could be another error in my implementation.
I suspect this as the build-in Neumann boundary works accordingly, but this basic fixed boundary does not, based on equivalent Matlab code.

The problem is described by Fick's second law:
$\frac{\partial u\left(x,t\right)}{\partial t}=D\frac{\partial u\left(x,t\right)}{\partial x^2}$

It is assumed that the concentration at the boundary is kept constant
$u_L = 0$ and $u_R = 1$

import numpy as np
import matplotlib.pyplot as plt
from pde import CartesianGrid, ScalarField, PDEBase, FieldCollection, MemoryStorage, PlotTracker, PDE

# Define the spatial and temporal domain
x_min, x_max = 0, 1
nx = 100
t_max = 0.5
nt = 20

# Create a Cartesian grid
grid = CartesianGrid([[x_min, x_max]], [nx])

# Define the initial condition
initial_value = 1
field = ScalarField(grid, initial_value)

# Define the boundary conditions of u
bc_left = {"value": "0"}
bc_right = {"value": "1"}
bc = [bc_left, bc_right]

# Applying boundary conditions and Laplacian
field.laplace(bc)

# Diffusion constant
diffusivity = 1

# Specifying the PDE for diffusion
eq = PDE({"u": f"{diffusivity}*laplace(u)"}, bc=bc)

# Create storage to save the solution
storage = MemoryStorage()

# Solve the PDE
result = eq.solve(field, t_range=t_max, dt=t_max/nt, tracker=storage.tracker(interval=0.001))

# Extract data from storage for plotting
times = np.linspace(0, t_max, len(storage.data))
data = np.array(storage.data)

# Convert the stored data into a format suitable for plotting
solution = data.reshape((len(storage.data), nx))

# Create a 2D meshgrid for plotting
X, T = np.meshgrid(grid.axes_coords[0], times)

# Plot the solution
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, T, solution, cmap='viridis')
fig.colorbar(surf)
ax.set_xlabel('Distance x')
ax.set_ylabel('Time t')
ax.set_zlabel('Solution u')
plt.title('Diffusion on a line')
plt.show()

The output plot is pictured below:

I have previously solved this problem in Matlab, and here the output looks like this:

I would very much appriacte your help on what could be causing this difference! Thanks!

[david-zwicker]

Your code runs into the classical von Neumann instability. To avoid this, it is better to use an adaptive time stepper by replacing the central line that solves the PDE by result = eq.solve(field, t_range=t_max, adaptive=True, tracker=storage.tracker(interval=0.001))

[GijsBerkenbosch]

Thank you very much for your swift reply David, I was not aware of this instability!

[GijsBerkenbosch]

Hi @david-zwicker,

I am still working on replicating Matlab code using the py-pde package, but have encountered another problem.

This problem focusses on the Schnakenberg model which involves 2 chemical species X and Y which undergo the following reactions:

The 2 PDE of the reaction diffusion system are given by:

We now look for the concentration of X and Y along a line which we can model as an uni-dimensional problem, stretching between $a$ and $b$ (thus $a\leq x \leq b$).
Since we have two coupled PDE we have two initial conditions, one for $u(x,t)$ and one for $v(x,t)$. The following initial conditions will be applied that generate, hopefully, nice spatial patterns
$u(x,0) = 0$
$v(x,0) = 0$

We take the system to be isolated: no molecules enter of exit the system. We then need Neumann's conditions applied at both ends of the line.
$\frac{\partial}{\partial x} v(x,t) = 0$
$\frac{\partial}{\partial x} u(x,t) = 0$ on $x=a,b$

`import numpy as np
import matplotlib.pyplot as plt
from pde import CartesianGrid, ScalarField, PDE, MemoryStorage, FieldCollection

Define the spatial and temporal domain

x_min, x_max = 0, 1
nx = 500
t_min, t_max =0, 10

Create a Cartesian grid

grid = CartesianGrid([[x_min, x_max]], [nx])

Define initial conditions

u0 = ScalarField(grid, 0)
v0 = ScalarField(grid, 0)
state = FieldCollection([u0, v0])

Define the boundary conditions

boundary_conditions = "neumann"

Model parameters

a = 0.2
b = 0.5
g = 70
d = 100

Define the PDEs

eq = PDE({
"u": f"laplace(u) + {g} * ({a} - u + u2 * v)",
"v": f"{d} * laplace(v) + {g} * ({b} - u2 * v)"
}, bc=boundary_conditions)

Create storage to save the solution

storage = MemoryStorage()

Solve the PDE

result = eq.solve(state, t_range=t_max, adaptive=True, tracker=storage.tracker(interval=0.1))

Extract data from storage for plotting

times = np.linspace(0, t_max, len(storage))
data = np.array(storage.data)
data_u = data[:, 0, :]
data_v = data[:, 1, :]

Convert the stored data into a format suitable for plotting

solution_u = data_u.reshape((len(storage), nx))
solution_v = data_v.reshape((len(storage), nx))

Create a 2D meshgrid for plotting

X, T = np.meshgrid(grid.axes_coords[0], times)

Plot the solution

fig, axs = plt.subplots(2, 1, figsize=(10, 8))

Plot u

cax1 = axs[0].imshow(solution_u, aspect='auto', origin='lower', cmap='jet', extent=[x_min, x_max, t_min, t_max])
axs[0].set_ylabel('time')
axs[0].set_title('Example C: Schnakenberg-Touring patterns')
fig.colorbar(cax1, ax=axs[0])

Plot v

cax2 = axs[1].imshow(solution_v, aspect='auto', origin='lower', cmap='jet', extent=[x_min, x_max, t_min, t_max])
axs[1].set_xlabel('space x')
axs[1].set_ylabel('time')
fig.colorbar(cax2, ax=axs[1])

plt.tight_layout()
plt.show()
`
The obtained figure:

While the matlab code looks like this:

Both the solutions for u and v are different from Matlab. U appears constant in space x, while this is not the case in Matlab, and v has a slightly different pattern with higher values. Do you have any idea what could be causing this problem? Thank you very much!

[david-zwicker]

I'm afraid I don't have the bandwidth to investigate this in depth. You might be interested in the plot_kymograph function of py-pde, though.