Improve solution speed of PDE with nonlocal right-hand side

General
1

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] 

mesh = UnitSquareMesh.create(mesh_resolution,mesh_resolution, CellType.Type.quadrilateral)
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 * 0.1 * inner(grad(Y),grad(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-------")
2

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

3

Thank you very much. I was able to shorten the runtime from 14-15s to 0.9-1s per time step. Amazing, thank you :slight_smile: 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] 

mesh = UnitSquareMesh.create(mesh_resolution,mesh_resolution, CellType.Type.quadrilateral)
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 * 0.1 * inner(grad(Y),grad(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-------")