Time dependent Neumann BC leads to Calling FFC just-in-time (JIT) compiler at every timestep

1

Hi all,

I’m trying to implement a time-dependent Neumann BC. For instance, at every time step,

L_final = L + L_BC

where L is the same for every timestep and L_BC = f(t). My current challenge is because of this L_BC, at every timestep, whenever I call

solve(a == L_final, u, bc)

It will trigger “Calling FFC just-in-time (JIT) compiler” which is super time-consuming. I’m running for 100k time step.

I’m wondering if anyone has experience in this kind of implementation that avoids recompiling at every time step?

Thanks.

2

You could do something along the lines of

L_BC = Expression("t", t=0.0, degree=0) # either pass the element or the degree

# create the DirichletBC object with L_BC... and the list of bcs `bc`

# assuming T is your array of time-steps, say `numpy.linspace(0.,1.e6, 100000)`
for t in T:
    L_BC.t = t
    solve(....)

You could also equivalently use Constant and it’s corresponding assign method to do something similar, although its hard to tell without the problem and a minimal example.

1 Like
3

Thank you for your reply. I can be more specific on this.

Here is the simplified code I wrote.

L_BC = 0
for i in range(len(areas)):
      L_BC -= dt*qn[i]*v[0] * ds_p(i+1) 
      L_BC -= dt*qn[i]*v[1] * ds_c(i+1) 

L_final = L + L_BC
solve(a == L_final, u, bc)

Basically, qn is an array and changes for every time step. ds_p and ds_c are different subdomains of the Neumann boundaries.

Thank you for your help. Really appreciated.

4

Ah I see, I completely missed reading that it is a time dependent Neumann BC and not Dirichlet. What I wrote above is for Dirichlet BC. Nevermind !

As for your question, is the solve inside or outside the loop? Because with the code you have posted there should be nothing that requires any recompilation of the forms. Can you create a minimal working example instead? See Read before posting: How do I get my question answered?

5

I have created a simplified version of my code.

from fenics import *
import numpy as np

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)


def right(x, on_boundary):
    return near(x[0], 1)
def btm(x, on_boundary):
    return near(x[1], 0)
def top(x, on_boundary):
    return near(x[1], 1)


bc_right = DirichletBC(V, u_D, right)
bc_btm = DirichletBC(V, u_D, btm)
bc_top = DirichletBC(V, u_D, top)

bcs = [bc_right, bc_btm, bc_top]

left_bc_markers  = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
sub_domains = np.linspace(0, 1, 11)
for i in range(len(sub_domains)-1):
    class left_bc(SubDomain):
        def inside(self, x, on_boundary):
            return x[1] <= sub_domains[i+1] and x[1] >= sub_domains[i]  and near(x[0], 0)
    bx0 = left_bc()        
    bx0.mark(left_bc_markers,i+1)

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0


for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    ds_sub = Measure("ds", domain = mesh, subdomain_data=left_bc_markers)
    qn = np.random.default_rng().uniform(1,10,len(sub_domains)-1)
    print(qn)
    L_BC = 0
    for i in range(len(sub_domains)-1):
        L_BC += dt*v*qn[i]*ds_sub(i+1)

    L_Final = L + L_BC
    # Compute solution
    solve(a == L_Final, u, bcs)

    # Update previous solution
    u_n.assign(u)

The key idea is the Neumann boundary depends on the location discretely, so I divided the boundary into N subdomains with each subdomains having a value of qn.

6

How about wrapping it as a constant before entering the time-loop like I mentioned above? Does that work?

from fenics import *
import numpy as np

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)


def right(x, on_boundary):
    return near(x[0], 1)
def btm(x, on_boundary):
    return near(x[1], 0)
def top(x, on_boundary):
    return near(x[1], 1)


bc_right = DirichletBC(V, u_D, right)
bc_btm = DirichletBC(V, u_D, btm)
bc_top = DirichletBC(V, u_D, top)

bcs = [bc_right, bc_btm, bc_top]

left_bc_markers  = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
sub_domains = np.linspace(0, 1, 11)
for i in range(len(sub_domains)-1):
    class left_bc(SubDomain):
        def inside(self, x, on_boundary):
            return x[1] <= sub_domains[i+1] and x[1] >= sub_domains[i]  and near(x[0], 0)
    bx0 = left_bc()        
    bx0.mark(left_bc_markers,i+1)

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0

ds_sub = Measure("ds", domain = mesh, subdomain_data=left_bc_markers)
qn = np.random.default_rng().uniform(1,10,len(sub_domains)-1)
qn_wrapped = Constant(qn)
L_BC = sum(dt*v*qn_wrapped[i]*ds_sub(i+1) for i in range(len(sub_domains)-1))
# for i in range(len(sub_domains)-1):
#     L_BC += dt*v*qn[i]*ds_sub(i+1)

L_Final = L + L_BC

for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    
    qn_wrapped.assign(Constant(np.random.default_rng().uniform(1,10,len(sub_domains)-1)))
    # Compute solution
    solve(a == L_Final, u, bcs)

    # Update previous solution
    u_n.assign(u)
3 Likes
7

This works like a charm! Thanks!