1
$\begingroup$

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)
$\endgroup$
  • 1
    $\begingroup$ In newer versions of fenics, the parameter "reuse_factorization" seems to be removed. That could be causing the segmentation problem. Comment that line and try. $\endgroup$ – cpraveen Oct 9 at 8:44
  • $\begingroup$ Great. After the commenting, it works: much faster, less memory usage. Thanks so much! $\endgroup$ – duccio Oct 9 at 8:59
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I tried as suggested, but I obtain a segmentation fault. I'll edit my original answer adding a complete replicable short code. $\endgroup$ – duccio Oct 9 at 8:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.