Execution time with time-dependent sources

dolfinx
1

Hello,

I’m working on a transient thermo-mechanical problem and I’ve noticed that the execution time increases significativally when I substitute the constant source term with a time dependent source term.

With the constant source (declared as f = fem.Constant(msh, ScalarType((0.0, 0.0))) ), execution takes around 1.5s, whereas with the time dependent source, it takes around 117.6s. Therefore, I’m wondering if there is a more efficient way of defining or updating this function.

The strategy I’m using is to define the source term as a function f_exp(x,t) and using interpolate for uptading at each step t, ie.,

f.interpolate(partial(f_exp_1, t)).

Indeed, I’ve noticed that it is the interpolation itself that is consuming a lot of time. As illustrated on the MWE bellow, interpolationg for a fixed time step takes around 6 seconds, even with a small mesh (20x20). Hence, when executing the transient problem, with several time steps, the execution time is very high.

I wonder if anyone could help me understand why this interpolation is taking so long. I suspect it has something to do with the way I define this long expression, since the second one (simpler) takes substantially less time.

Exists on my machine were:

Execution time 1: 6.352968215942383
Execution time 2: 0.002003192901611328

Any suggestions? I appreciate very much any help.

# Auxiliary packages
from petsc4py import PETSc # Linear algebra backend
from mpi4py import MPI     # MPI 
import numpy as np
from matplotlib import pyplot as plt

# Import dolfinx and ufl modules
from dolfinx import geometry
from dolfinx import fem
from dolfinx import io
from dolfinx import la
from dolfinx.mesh import (CellType, create_unit_square, locate_entities, exterior_facet_indices, create_rectangle, meshtags)
from ufl import sin, SpatialCoordinate, FiniteElement, VectorElement, MixedElement, TestFunction, TrialFunction, split, Identity, Measure, dx, ds, grad, nabla_grad, div, dot, inner, tr, as_vector, FacetNormal
from petsc4py.PETSc import ScalarType
from functools import partial 

# Import dolfinx modules and ufl
from dolfinx.io import gmshio
import dolfinx.cpp as _cpp
import time

Lx = 1.0; Ly = 1.0; N = 20

msh = create_rectangle(comm=MPI.COMM_WORLD,
                    points=((0.0, 0.0), (Lx, Ly)), n=(N, N),
                    cell_type=CellType.triangle, diagonal = dolfinx.cpp.mesh.DiagonalType.crossed)

E, nu = 30e9, 0.25                                                       # Young modulus, Poisson ratio
mu    = fem.Constant(msh, ScalarType(E/(2.0*(1.0 + nu))))                # Lamé constants
lbd  = fem.Constant(msh, ScalarType(E*nu/((1.0 + nu)*(1.0 - 2.0*nu))))
alpha    = fem.Constant(msh, ScalarType(12e-6))                          # thermal expansion     

### Function spaces ###
degree_u = 1; degree_T = 1
el_u = VectorElement("CG", msh.ufl_cell(), degree_u, dim=2) 
el_T = FiniteElement("CG",  msh.ufl_cell(), degree_T)
el_m  = MixedElement([el_u , el_T]) 
W = fem.FunctionSpace(msh, el_m) 
U, _ = W.sub(0).collapse() 
V, _ = W.sub(1).collapse() 

def f_exp_1(t, x):
    return (np.full_like(x[0], (np.pi**2)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1])- 
                         (np.pi**2)*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])-alpha*np.pi*np.exp(-np.pi*t)*np.cos(np.pi*x[0])), 
            np.full_like(x[0], (np.pi**2)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1])- 
                         (np.pi**2)*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])-alpha*np.pi*np.exp(-np.pi*t)*np.sin(np.pi*x[1])))

ti = time.time()
f = fem.Function(U)
f.interpolate(partial(f_exp_1,0)) 
tf = time.time()

print('Execution time 1: '+str(tf - ti))


def f_exp_2(t, x):
    return (np.full_like(x[0], np.pi*np.exp(-np.pi*t)*np.cos(np.pi*x[0])), 
            np.full_like(x[0], np.pi*np.exp(-np.pi*t)*np.sin(np.pi*x[1])))

ti = time.time()
f = fem.Function(U)
f.interpolate(partial(f_exp_2,0)) 
tf = time.time()

print('Execution time 2: '+str(tf - ti))```
2

I get

Execution time 1: 0.9482619762420654
Execution time 2: 0.0005140304565429688

for your code. Your f_exp_1 function is a lot more complicated than f_exp_2 if I put timers inside, e.g.

def f_exp_1(t, x):
    ti = time.time()
    y = (np.full_like(x[0], . . . . 
    tf = time.time()
    print('Execution time F1: '+str(tf - ti))
    return y

I get

Execution time F1: 0.9475531578063965
Execution time F2: 8.0108642578125e-05

indicating that nearly all the time is in the NumPy calls. You could try to re-arrange the NumPy calls, or maybe Numba can help.

If your time-dependent source has the form f = g(t) h(x), you could interpolate h at the start, and update the degrees-of-freedom of f directly at each time step, e.g.

h = fem.Function(U)
h.interpolate(f_exp)

f = fem.Function(U)
for i in range(20):
    t = i * dt
    f.x.array[:] =  t * h.x.array
1 Like
3

This does not make sense to me, as you shouldn’t need to use full_like at all. If you use full_like on an array, you get a matrix back, see: numpy.full_like — NumPy v1.26 Manual

Tou should just keep

def f_exp_1(t, x):
    return ((np.pi**2)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1])- 
                         (np.pi**2)*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])-alpha*np.pi*np.exp(-np.pi*t)*np.cos(np.pi*x[0]))
, 
 (np.pi**2)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1])- 
                         (np.pi**2)*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])-alpha*np.pi*np.exp(-np.pi*t)*np.sin(np.pi*x[1]))
4

Thank you @garth and @dokken for you quick response and suggestions. It was very helpfull.

Still regarding the time update strategy, I’ve seen in some tutorials the use of classes to define time dependent boundary conditions (for exemple in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial, for the inlet velocity). Is it a good practice to always define classes for time depedent functions?

5

I thought your use of partial was neat. I don’t see any problem.

6

I would say it depends on what you need to store for your expression. As you only store the time, using partial is fine.