solve multiple equations at once

Question

I am looking for an efficient way to solve multiple equations at once from python. I want to solve

$$ \frac{dI(x,E)}{dx} + \alpha(E)I(x,E) = 0 $$

For a given E, I can solve it easily like with the code snippet below. My problem is that if I want to solve it for several E's I need a for loop which is very inefficient in python. The solution I am after can be written as

$$ I(x) = \int{ I(x,E)dE } $$

Any idea would be appreciated.

from dolfin import *
m1 = IntervalMesh(10,0.0, 1.0)
Q1 = FunctionSpace(m1,'CG',1)
v1 = TestFunction(Q1)
du1 = TrialFunction(Q1)
dx = Measure('dx')

def inside_left_m1(x, on_boundary):
    return near(x[0],0.0)
left_m1_domain = AutoSubDomain(inside_function=inside_left_m1)
left_m1_bc = DirichletBC(Q1,1.0,left_m1_domain)
alpha = Constant(1.0)

F = du1.dx(0)*v1 + alpha(1.0)*du1*v1
F *= dx

a = lhs(F)
L = rhs(F)

u1 = Function(Q1)
parameters = LinearVariationalSolver.default_parameters()
parameters['linear_solver'] = 'lu'

solve(a == L, u1, [left_m1_bc,],
        solver_parameters=parameters)

plot(u1)

Answer 1

If I understand you right, you're able to formulate your problem as a variational problem involving double integral for unknown function of two variables. So this is a 2D problem but certainly you should exploit a structure of your problem to split it into two one-dimensional problems for efficiency. I have a suspicion that you're going into premature optimization with eliminating python loops but what you need to do is to eliminate a need for calling JIT compiler for every value of E.

Comments

Sorry I think E should be seen more as a parameter unless there is a way to treat it as a variable in dolfin. Here is a working example of what I am trying to achieve. The problem is that as I increase the number of energies,intensities will not scale well in speed (replace 101 with 2000). Maybe there is a way to speed this up using some type of MixedFunctionSpace formulation? Otherwise I would have to use c++ but my code is in python.

from dolfin import *
import numpy as np
import pylab as pb

m1 = IntervalMesh(10,0.0, 1.0)
Q1 = FunctionSpace(m1,'CG',1)

energies = np.linspace(1.0,10,101)
intensities = np.linspace(0.1, 1.0,101)

v1 = TestFunction(Q1)
du1 = TrialFunction(Q1)

dx = Measure('dx')

def inside_left_m1(x, on_boundary):
    return near(x[0],0.0)

left_m1_domain = AutoSubDomain(inside_function=inside_left_m1)

alpha = Constant(1.0)
light_intensity = {}

for energy, intensity in zip(energies,intensities):
    F = du1.dx(0)*v1 + alpha(energy)*du1*v1
    F *= dx

    a = lhs(F)
    L = rhs(F)

    u1 = Function(Q1)

    parameters = LinearVariationalSolver.default_parameters()

    parameters['linear_solver'] = 'lu'

    left_m1_bc = DirichletBC(Q1,intensity,left_m1_domain)

    solve(a == L, u1, [left_m1_bc,],
        solver_parameters=parameters)

    light_intensity[energy] = u1

---

Chaffra, I didn't try the code but I can see the problem. You use alpha(energy) which is float in your form. (Constant.__call__ returns float!). Hence your form is recompiled every iteration. Try Constant(alpha(energy)) or better move your form out of loop and just use Constant.assign() method to change its parameters every iteration.

For optimization you could also assemble system matrix and system vector only once. Then you just create copies of matrix and vector, multiply vector by your alpha(energy) and apply bcs every iteration.