Hello everyone,

I’d like to use FEniCS to calculate an integral over a (discretized) domain. The integral in question is the weak form of a PDE (e.g. a diffusion PDE for now), but I do not want to solve a finite element problem. I, instead, want to plug in (discretized) functions for what is usually a Test- and TrialFunction and calculate a single scalar value (i.e. residual) using the weak form. To reemphasize, I do not need a solver, simply the FEniCS function “assembly” to evaluate the integrals in the weak form. Also, I need the derivatives of the weak form w.r.t. the function I plug in.

Here is my problem: I need to evaluate this integral in the order of 10^5 - 10^6ish times, and this evaluation is the bottle neck of my project. I guess my implementation is suboptimal and I have no experience what so ever with speeding up things in FEniCS. I have a minimal working example attached.

Please consider the following two cases:

- The form of my integral remains the same, but I need to plug in new numbers for all my (discretized) input functions,
- The form of the weak form remains the same and some of the input functions remain constant, while I need to plug in new numbers for other input functions.

Is there any thing that can be done to speed up this process?

Thanks in advance for your help!

Best, Vincent

Here is the minimal working example:

```
import fenics as df
import numpy as np
#%% Setup
# Basic setup, where there may be multiple function spaces
mesh = df.UnitSquareMesh(8, 8)
V = df.FunctionSpace(mesh, 'CG', 1)
Vc = df.FunctionSpace(mesh, 'DG', 0)
# create some functions for my formula
u = df.Function(V)
v = df.Function(V)
k = df.Function(Vc) # <-- (potentially) different function space
# initialize some random values
v.vector().set_local(np.random.rand(V.dim()))
k.vector().set_local(np.random.rand(Vc.dim()))
# some arbitrary load
x = df.SpatialCoordinate(mesh)
f_load = 1 * x[0] + 2 * x[1]
#%% Part I have to run multiple times
some_large_number = 50 # <-- This is usually about 10^5
for i in range(some_large_number):
# calculate new values for my functions (here: just random)
# case a) some functions remain constant (e.g. v, k), and some change (e.g. u), as shown in this minimal example
# case b) all functions (u, v, k) change
u.vector().set_local(np.random.rand(V.dim()))
# evaluate my formula res(u, v, k), which returns a scalar (!) residual
# >>> first call: 3e-4 s, successive calls: 2e-4 s
res_form = k * df.inner(df.grad(u), df.grad(v)) * df.dx - f_load * v * df.dx
# >>> first call: 1e-3 s, successive calls: 4e-4 s
res = df.assemble(res_form)
# calculate the respective derivatives (here: just u for clarity)
# >>> first call: 1e-3 s, successive calls: 7e-4 s
dres__du = np.asarray(df.assemble(df.derivative(res_form, u)))
```