Stopping Just-In-Time Compiler every loop of Monte Carlo Simulation

Question

Sure this question highlights my lack of knowledge of about how the compilers work.

We are trialling out using Fenics for using solve Monte Carlo problems, solving Poisson equations with random coefficients.

The problem is in loops where a new random field is define each time, when you run it the computational cost is killed by the 'Just-In-Time' Compiler.

Here would be a simple example, where I use the poisson example and randomise the right hand side

Any help would be great.

Thank you

Tim Dodwell, University of Exeter.

from dolfin import *

import random
import numpy as np

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
u0 = Function(V)
bc = DirichletBC(V, u0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
               degree=1)
g = Expression("sin(5*x[0])", degree=1)
a = inner(grad(u), grad(v))*dx()

for  i in range(0,2):

    L = random.random()*f*v*dx() + g*v*ds()

    # Define function for the solution
    u = Function(V)

    # Define goal functional (quantity of interest)
    M = u*u*dx()

    # Define error tolerance
    tol = 1.e-5

    # Solve equation a = L with respect to u and the given boundary
    # conditions, such that the estimated error (measured in M) is less
    # than tol
    problem = LinearVariationalProblem(a, L, u, bc)
    solver = AdaptiveLinearVariationalSolver(problem, M)
    solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
    solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
    solver.solve(tol)

    solver.summary()

Magne Nordaas · Answer 1 · Oct 25, 2016

You can avoid the excessive jitting using Constant scalars. For example,

r = Constant(random.random())
L = r * f * v * dx() + g * v * ds()

answered Oct 25, 2016 by Magne Nordaas FEniCS Expert (13,820 points)

jenom · Answer 2 · Nov 14, 2016

As Magne Nordaas pointed out in another answer you have to avoid the recompilation. One method is to use a Constant() as he pointed out. If the random part is a field however this will not work. You have two options:

Use a precompiled Expression with parameters that change. Use either the C++-Code interface (Variant 2 in Expression) for the whole class or the parameters (Variant 1 in Expression), that you can extent according to the interface. Avoid Variant 3 in Expression as it will be much slower.
Interpolate your random field onto the common mesh with some fixed element and degree before including it into the form.

It is also possible to compute the samples in parallel with FEniCS though it is not best practice and needs some extra care. Let me know if you want to know more.

Best,
Johannes

Stopping Just-In-Time Compiler every loop of Monte Carlo Simulation

3 Answers