Dealing with PDEs on an Infinite Domain

General
1

Hi everyone,

I’m using FEniCS 2019.1.0 to solve the following system of differential equations

\frac{\partial c}{\partial t} = \nabla^2c - c + \delta(\mathbf{x}; \mathbf{x}_a(t) ) \tag{1}
\frac{d\phi}{dt} = \Lambda \nabla c \cdot \begin{pmatrix} -\sin \, \phi(t) \\ \cos\,\phi(t) \end{pmatrix} \tag{2}
\frac{d\mathbf{x}_a}{dt} = v\begin{pmatrix} \cos\,\phi(t) \\ \sin\,\phi(t) \end{pmatrix} \tag{3}

where c represents the concentration of a chemical substance and \mathbf{x}_a(t) is the position of a particle at time t. In this case, (1) has an analytical solution,

c(\mathbf{x},t) = \frac{1}{4\pi} \int_0^{t - \epsilon t} \frac{e^{-(t - \tau)} }{t - \tau} \exp\biggl( - \frac{ [ \mathbf{x} - \mathbf{x}_a(t) ]^2} {4(t - \tau)} \biggr) d\tau \tag{4}

where \epsilon t is a small interval used to regularize the above.

Unfortunately, the numerical results I’m getting don’t coincide with the analytical solution given by (4). Solving (1) - (3) on the unit square, with \mathbf{x}_a(0) = (0.5, 0.5)\; \phi(0) = \pi/3, using a time-step \Delta t = 0.0005 with mesh-size h = 1/50, after 1000 time-steps I obtain the following numerical result:

num_soln

Conversely, with the analytical solution, after 1000 time-steps (using the same value of \Delta t), I have the following -

anal_soln

One salient difference between the two programs (i.e. the one using the analytical solution and the program using the FEM solution) is in how c decays spatially. Indeed, in the case of the former, after 2 time-steps c decays by about 6 orders of magnitude going a distance 10h away from the point source at \mathbf{x}_a. Conversely, in the case of the numerical solution – again, after 2 time-steps – if one looks a distance 10h away from the point source, c decays by only about 3 orders of magnitude.

Based on what I’ve seen thus far, making the mesh finer and reducing the time-step have no effect on the accuracy of the numerical solution. I can think of two possible causes for the pathological behavior of the numerical solution:

  1. Boundary conditions. Indeed, (1) obeys boundary conditions at \infty, whereas I am forced to solve the differential equations on a finite domain. I have experimented with both homogeneous Neumann and Dirichlet boundary conditions and neither leads to accurate results. One would perhaps intuit that making the computational domain larger would reduce the effect of any imposed boundary conditions and completely negate any such effects in the limit as the domain becomes arbitrarily large. Again, this does not seem to be the case. When I make the domain larger, c still decays far too slow and the results I get don’t coincide with the analytical solution.

  2. Could this possibly be a fundamental limitation of the PointSource method?

Can anyone think of another possible cause of this bad behavior? Does anyone have any suggestions for remedies? Find attached (first) the program I’m using to compute the FEM solution and below it the program I’m using to compute the analytical solution.

from fenics import *
import numpy as np
import sys
import time

start_time = time.time()

Lambda = 4000

dt = 5e-4 # time-step
Nt = 1000

a, b, c, d = 0, 1, 0, 1
# Create mesh 
nx = ny = 50
#mesh = RectangleMesh(Point(a, b), Point(c, d), nx, ny)
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Used to plot 'c'
x_vec = np.linspace(a, b, nx + 1)
y_vec = np.linspace(c, d, ny + 1)
[X, Y] = np.meshgrid(x_vec, y_vec)
C_fem = np.zeros( (nx+1, nx+1, Nt) )


c_0 = Expression('0',degree = 2)
c_D = Expression('0*x[0] + 0*x[1]', degree = 1)

def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, c_D, boundary)

c_n = interpolate(c_0, V)


c = TrialFunction(V)
v = TestFunction(V)

dt_ufl = Constant(dt)

a = (1 + dt_ufl)*c*v*dx + dt_ufl*dot( grad(c), grad(v) )*dx
L = c_n*v*dx
c = Function(V)

particlePosition = np.zeros( (Nt,2) )
particlePosition[0,0] = 0.5
particlePosition[0,1] = 0.5

phi = np.zeros(Nt+1)
phi[0] = np.pi/3
speed = 0.05

assembler = SystemAssembler(a, L, bc)
solver = LUSolver("mumps")
A = Matrix()
solver.set_operator(A)
assembler.assemble(A)


class Projector():
    def __init__(self, V):
        u = TrialFunction(V)
        v = TestFunction(V)
        a = inner(u, v) * dx
        self.A = assemble(a)
        self.u = Function(V)
        self.solver = LUSolver(self.A)

    def __call__(self, f):
        v = TestFunction(self.u.function_space())
        b = assemble(inner(f, v)*dx)
        self.solver.solve(self.u.vector(), b)
        return self.u


V_vec = VectorFunctionSpace(mesh,"CG",1)

t = 0
for n in range(Nt-1):
    # Update time
    t += dt

    # Point where I am evaluating the point source. If n \neq 0, I evaluate
    # the point source at x_a(n-1) - the previous position of the particle
    # to regularize the integral in equation (1)
    if n == 0:
        pos_curr = particlePosition[n,:]
    else:
        pos_curr = particlePosition[n - 1,:]
    pt = Point(pos_curr[0],pos_curr[1])


    # Compute solution
    b = assemble(L)
    delta = PointSource(V, pt, dt)
    delta.apply(b)
    # bc.apply(b)
    solve(A, c.vector(), b)

    c_n.assign(c)
    
    # for i in range(len(x_vec)):
    #     for j in range(len(y_vec)):
    #         C_fem[i, j, n] = c( (x_vec[i], y_vec[j]) )

    projector = Projector(V_vec)
    gradc = projector(grad(c))
    grad_c = gradc(( pos_curr[0],pos_curr[1] ))
    dphi_dt = Lambda*(-grad_c[0]*np.sin(phi[n])
                 + grad_c[1]*np.cos(phi[n]))
    phi[n+1] = dt*dphi_dt + phi[n]
        

    if n == 0:
        particlePosition[n+1,0] = speed*np.cos(phi[n])*dt + particlePosition[n,0]
        particlePosition[n+1,1] = speed*np.sin(phi[n])*dt + particlePosition[n,1]
    else:
        particlePosition[n+1,0] = speed*np.cos(phi[n+1])*dt + particlePosition[n,0]
        particlePosition[n+1,1] = speed*np.sin(phi[n+1])*dt + particlePosition[n,1]  

#%% Plot the trajectory of the particle
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
[fig,ax] = plt.subplots()
ax.plot(particlePosition[:,0],particlePosition[:,1])
ax.set_aspect('equal')
ax.set_xlabel(r"$x$",fontsize=20)
ax.set_ylabel(r"$y$",fontsize=20)
ax.text(0.504, 0.500, r"$\Lambda = %g$" %Lambda, color="k", fontsize=20)
ax.text(0.508, 0.500, r"$v = %g$" %speed, color="k", fontsize=20)

from fenics import * # analysis:ignore
import numpy as np
import matplotlib.pyplot as plt
import time

Lambda = 4000

dt = 5e-4 # Time-step
Nt = 1000

a, b, c, d = 0, 1, 0, 1
nx = ny = 50

# For plotting c
x_vec = np.linspace(a, b, nx + 1)
y_vec = np.linspace(c, d, ny + 1)
[X, Y] = np.meshgrid(x_vec, y_vec)
C_analytical = np.zeros( (nx+1, ny+1, Nt) )

def get_chem_field(r_a, x0, y0, n, dt):
    integ = 0
    
    r = np.array([x0, y0])
    for i in range(n):
        
        dist = r - r_a[i]
        dist_squared = np.linalg.norm( dist )**2

        delta_t = n*dt - i*dt
        integ += (1/delta_t)*np.exp(-delta_t)*np.exp(-dist_squared/(4*delta_t))
    c_field = (1/4*np.pi)*integ*dt
    return c_field

def get_chem_grad(r_a, n, dt):
    h = 1e-4
    x_a, y_a= r_a[n,0], r_a[n,1]
    c_x0_plus_h = get_chem_field(r_a, x_a + h, y_a, n, dt)
    c_x0_minus_h = get_chem_field(r_a, x_a-h, y_a, n, dt)
    
    c_x = (c_x0_plus_h - c_x0_minus_h)/(2*h)
    
    c_y0_plus_h = get_chem_field(r_a, x_a, y_a + h, n, dt)
    c_y0_minus_h = get_chem_field(r_a, x_a, y_a - h, n, dt)
    
    c_y = (c_y0_plus_h - c_y0_minus_h)/(2*h)
    
    
    return np.array( [c_x, c_y] )

# Create arrays for particle position and angular orientation.
# Populate these arrays with initial values.
particlePosition = np.zeros((Nt,2))
particlePosition[0,0] = 0.5
particlePosition[0,1] = 0.5

phi = np.zeros(Nt+1)
phi[0] = np.pi/3
speed = 0.05

# Time-stepping
t = 0
grad_c = np.zeros((Nt,2))

for n in range(Nt-1):

    # Compute solution  
    
    grad_c[n,:] = get_chem_grad(particlePosition, n, dt)
    dphi_dt = Lambda*(-grad_c[n,0]*np.sin(phi[n])
                 + grad_c[n,1]*np.cos(phi[n]))
    phi[n+1] = dt*dphi_dt + phi[n]
    
    for i in range(len(x_vec)):
        for j in range(len(y_vec)):
            C_analytical[i, j, n] = get_chem_field(particlePosition,
                                        x_vec[i], y_vec[j], n, dt)
        
    if n == 0:
        particlePosition[n+1,0] = speed*np.cos(phi[n])*dt + particlePosition[n,0]
        particlePosition[n+1,1] = speed*np.sin(phi[n])*dt + particlePosition[n,1]
    else:
        particlePosition[n+1,0] = speed*np.cos(phi[n+1])*dt + particlePosition[n,0]
        particlePosition[n+1,1] = speed*np.sin(phi[n+1])*dt + particlePosition[n,1]
        

plt.rcParams['text.usetex'] = True
[fig,ax] = plt.subplots()
ax.plot(particlePosition[:,0],particlePosition[:,1], color="k")
ax.set_aspect('equal')
ax.set_xlabel(r'$x$',fontsize=20)
ax.set_ylabel(r'$y$',fontsize=20)

ax.text(0.501, 0.5005, r"$\Lambda = 4000$", color="k", fontsize=15)
ax.text(0.503, 0.5005, r"$v = 0.05$", color="k", fontsize=15)

end = time.time()
print('modified, time = ', end-start)
2

Without being intimately familiar with your specific problem, I have a couple of thoughts:

  • It looks like you’re using an explicit first-order time stepping approach? Try changing it to an implicit one (or alternatively try a higher-order explicit time stepping scheme). I say this because you first solve a linear system for c with an explicit first-order time stepping scheme, then you use the same scheme to compute \phi from c and then you again use an explicit time stepping scheme to compute \bf{x}_a from \phi. Those first-order errors all add up, so to get accurate answers you’ll likely need to use some higher-order time stepping schemes. Note that (2) is a nonlinear equation, so you’ll probably need to use a Newton solver for it.
  • How are you measuring the “accuracy of the numerical solution”? Are you just looking at the particle trajectory, or also comparing the error in c and/or \phi?
  • You point out that the analytical solution has (homogeneous?) boundary conditions at \infty. I wouldn’t expect that imposing homogeneous boundary conditions on some finite domain would lead to the same solution, since you’re fundamentally solving a different problem. I’m not very familiar with the common approaches to deal with this, but one thing you could do for now is to pose the exact solution on the finite domain boundary as an inhomogeneous boundary condition.