Fenics: solving the same PDE multiple times

Question

I am new to Fenics and just started reading the tutorial Solving PDEs in Python. For simplicity, we can refer to simplest example, page 17 (the linear poisson equation), despite not necessary.

My problem is: given a fixed PDE, solve it multiple times with different parameters. (For instance, changing the constant f the mentioned example).

My 30 seconds attempt (notation as in the link above):

samples = 10

def solver(x, y):
   #x is the parameter to vary,
   #y the variable containing my solution
    f = Constant(x)
    L = f*v*dx
    u = Function(V)
    solve(a == L, u, bc)
   #write part of the solution in y

for i in range(samples):
    #set x, y as appropriate
    solver(x,y)

does the job, but when samples goes around 100.000 (magnitude that I need), it needs a lot of time and memory (often the process ends to be killed). My question is: how can I speed it up?

My idea is: since the PDE is always the same, there must be a way e.g. to store the stiffness matrix one time, and then use it every time, reducing "solver" to a linear algebra multiplication. On the other hand I do not know how to interface properly with fenics. Any suggestion is kindly appreciated.

Solution The answer given by cpraveen solved the issue. Here there is a raw but effective code, possibly useful for other readers.

from fenics import *
import numpy as np

## Commets omitted for simplicity: standard introduction
mesh = UnitSquareMesh(8,8)
V = FunctionSpace(mesh, 'P', 1)
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
        return on_boundary

bc = DirichletBC(V, u_D, boundary)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

print("1")

# New suggested trick
A = assemble(a)
bc.apply(A)
solver = LUSolver(A)
#solver.parameters['reuse_factorization'] = True # <- deprecated!

print("2")

# ...and define the solver accordingly
def my_solver(x):
    f = Constant(x[0])
    b = assemble(f*v*dx)
    bc.apply(b)
    L = f*v*dx
    u = Function(V)
    solver.solve(u.vector(), b)

# Start multiple resolutions
print("Go!")
samples = 10
for i in range(10):
        my_solver([i])
        print(i)

Answer 1

The following approach should work.

A = assemble(a)
bc.apply(A)
solver = LUSolver(A)
solver.parameters['reuse_factorization'] = True # maybe not needed

Now you can solve for any right hand side f

b = assemble(f*v*dx)
bc.apply(b)
u = Function(V)
solver.solve(u.vector(),b)

This avoids assembly and LU factorization in every solve. Only LU solve is required to be done.

Of course this works if the problem size is small enough that LU is viable.