Optimizing Code

General
1

Hi guys,

I’m writing a program in FEniCS 2019.1.0 to solve the following coupled reaction-diffusion system

c_t = D_c\nabla^2c -kc + h\delta(\mathbf{x} - \mathbf{x}_a(t) ) \\ \frac{d}{dt}\phi = \frac{\kappa}{\gamma_R}\begin{pmatrix} -\sin\, \phi \\ \cos \phi \end{pmatrix}\cdot \nabla c \\ \frac{d}{dt}\mathbf{x}_a = v\begin{pmatrix} \cos\,\phi \\ \sin \, \phi \end{pmatrix}.

where D_c, \; k, \; h, \; \kappa, \; \gamma_R and v are positive constants. My current program (see below) runs and generates the correct results - however, it’s quite slow.

The one promising thing I tried was to avoid assembling the entire linear system (matrix and rhs vector), but this didn’t work. Does anyone have suggestions?

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

k = 10
t_ref = 1/k
D_c = 1e-9
l_ref = np.sqrt(D_c/k)
kappa=1000 # 100
Lx = 1
Ly = 1
T = 10*t_ref
Lambda = 4000
h = 1e-5

gamma_R = (h*kappa*t_ref**2)/(Lambda*t_ref**3)

dt = t_ref/200 # Time-step
Nt = int(round(T/float(dt)))

# Create mesh and define function space
nx = ny = 50
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

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.005

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)

# Define variational problem
c = TrialFunction(V)
v = TestFunction(V)

a = (1 + k*dt)*c*v*dx + D_c*dt*dot( grad(c),grad(v) )*dx
L = c_n*v*dx

# Time-stepping
c = Function(V)
t = 0
grad_c = np.zeros((Nt,2))
mod_grad_c = np.zeros(Nt)
for n in range(Nt-1):
    if n == 0:
        pos_curr = particlePosition[n,:]
    else:
        pos_curr = particlePosition[n-1,:]
    pt = Point(pos_curr[0],pos_curr[1])


    # Compute solution
    [A,b] = assemble_system(a,L,bc)
    delta = PointSource(V,pt,h*dt)
    delta.apply(b)
    solve(A, c.vector(), b)

    c_n.assign(c)

    V_vec = VectorFunctionSpace(mesh,"CG",1)
    gradc = project(grad(c),V_vec)
    grad_c[n,:] = gradc(( pos_curr[0],pos_curr[1] ))
    mod_grad_c[n] = np.sqrt(grad_c[n,0]**2 + grad_c[n,1]**2)
    dphi_dt = kappa/(gamma_R)*(-grad_c[n,0]*np.sin(phi[n])
                 + grad_c[n,1]*np.cos(phi[n]))
    phi[n+1] = dt*dphi_dt + phi[n]
        

    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.figure()
plt.plot(particlePosition[:,0],particlePosition[:,1])
plt.gca().set_aspect('equal')
2

First of all change this so that you use a reusable projection, see for instance Introduction to FEniCS - Darcy — Simula Summer School in Computational Physiology

You should also be able to re-use the matrix A for all time-steps by using the system-assembler directly.

assembler = SystemAssembler(a, L, bcs)
solver = LUSolver("mumps")
A = Matrix()
solver.set_operator(A)
assembler.assemble(A)
# Assemble b and add point source here
# .....

# Define solution vector x and solve
x_func = Function(V)

solver.solve(x_func.vector(), b)
1 Like
3

Thank you! This did result in a speed-up of my program.

Question - can I assemble the right-hand side vector by simply calling

b = assemble(L)

?

This is what I’ve been doing and it seems to work - however, it seems like this formulation doesn’t apply boundary conditions.

4

You need to apply the boundary condition explicitly after assembly if you do it in this way.