# Improve solution speed of PDE with nonlocal right-hand side

Hey. I have a nonlocal PDE of the form

\partial_t y - \nu \Delta y + y + aY + b \cdot \nabla Y + \int_\Omega Y =0

with Neumann boundary. Each solver step takes 14-15 seconds right now, which seems too slow, please find a minimal running example below. The main issue lies in calling solve(A==0,Y) at each time step, I guess, and not specifying any solver. Moreover, I have to assemble in each time step the term \int_\Omega Y in the nonlocal term, I am not sure if there is a better way to do it. Because of this, I am also not sure how to write the problem in a class NonlinearVariationalProblem (as in this question) since the problem has to be updated in each time step. I am happy for any help.

from dolfin import *
from fenics import *
from mshr import *
import time

dt = 1e-3
number_time_steps = 20
final_time        = number_time_steps*dt
mesh_resolution   = 100

class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 1-2*x[0]*x[1]

x = SpatialCoordinate(mesh)

P1   = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
P1space = FunctionSpace(mesh, "Lagrange", 1)
dY     = TrialFunction(P1space)
test   = TestFunction(P1space)
Y      = Function(P1space)
Y0     = Function(P1space)
Y_initial = InitialConditions(degree=1)
Y  = interpolate(Y_initial, P1space)
Y0 = interpolate(Y_initial, P1space)

a_fun   = Expression('-1.5 + x[0] - abs(sin(6*t+x[0]))', t=0, degree=2)
b1_fun  = Expression("x[0]+x[1]", degree=2)
b2_fun  = Expression("abs(cos(6*t) * x[0] * x[1])", t=0, degree=2)

weak = (Y * test * dx
- Y0 * test * dx
+ dt * Y * test * dx
+ dt * a_fun * Y * test * dx  #a
+ dt * b1_fun * Y.dx(0) * test * dx #b1
+ dt * b2_fun * Y.dx(1) * test * dx #b2
)
testdx = test * dx
Ydx    = Y * dx

t = 0.0

while (t < final_time):
t += dt
tt = time.time()
a_fun.t = t
b1_fun.t = t
b2_fun.t = t
assemble_Ydx = assemble(Ydx)
nonlocalterm = dt * assemble_Ydx * testdx
solve(weak + nonlocalterm  == 0, Y)
Y0.assign(Y)
elapsed = time.time() - tt
print("-------It took "+str(elapsed)+" seconds-------")


First of all, you should pre-allocate matrices and vectors for all terms (if not using nonlinearvariationalproblem.

Another thing that would make it possible to use nonlinearvariational problem is create a placeholder constant for assemble(Ydx) to avoid recompilation at every time step.

Please note that your problem is actually linearized with the assemble_Ydx term, and that you could use a linear solver

Thank you very much. I was able to shorten the runtime from 14-15s to 0.9-1s per time step. Amazing, thank you I am not sure if GMRES+ILU are the best solver parameters, but here I could try further stuff. For completion, my code:

from dolfin import *
from fenics import *
from mshr import *
import time

dt = 1e-3
number_time_steps = 20
final_time        = number_time_steps*dt
mesh_resolution   = 100

class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 1-2*x[0]*x[1]

x = SpatialCoordinate(mesh)

P1   = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
P1space = FunctionSpace(mesh, "Lagrange", 1)
Y     = TrialFunction(P1space)
test = TestFunction(P1space)
Y_n  = Function(P1space)
Y_   = Function(P1space)
Y_initial = InitialConditions(degree=1)
Y_n = interpolate(Y_initial, P1space)

a_fun   = Expression('-1.5 + x[0] - abs(sin(6*t+x[0]))', t=0, degree=2)
b1_fun  = Expression("x[0]+x[1]", degree=2)
b2_fun  = Expression("abs(cos(6*t) * x[0] * x[1])", t=0, degree=2)

a1 = (Y * test * dx
+ dt * Y * test * dx
)
Abf1 = assemble(a1)
b   = test * dx
Lbf2 = assemble(b)
Ydx = Y * dx

t = 0.0

while (t < final_time):
t += dt
tt = time.time()
a_fun.t = t
b1_fun.t = t
b2_fun.t = t
a2 = (dt * a_fun * Y * test * dx  #a
+ dt * b1_fun * Y.dx(0) * test * dx #b1
+ dt * b2_fun * Y.dx(1) * test * dx #b2
)
Ydx_assemble = assemble(Ydx)
Lbf1 = assemble(Y_n * test * dx)
Lbf  = Lbf1 + Ydx_assemble * Lbf2
Abf = Abf1 + assemble(a2)
solve(Abf,Y_.vector(),Lbf,"gmres", "ilu")
Y_n.assign(Y_)
elapsed = time.time() - tt
print("-------It took "+str(elapsed)+" seconds-------")