Product of two finite UFL functions is non-finite

General
1

Hello everyone,

I have been working on implementing k-\epsilon (k-epsilon) turbulence model in FEniCS (I have already created one post about it here, here is the link if anyone is interested: Implementing k-epsilon model in FEniCS - how to implement terms in variational formulation that depend on a solution? - variational formulation - FEniCS Project).

The question I have is not related to my previous post or to turbulence modeling in general, however I encountered the error during my work on the model. That is also the reason why my minimal working example is not so minimal, since I haven’t been able to recreate the error outside of the environment I am working in.

Now, to my problem: I have two UFL functions f_1(k, \epsilon) and f_2(k, \epsilon). Here k and \epsilon are FEniCS Functions (computed solutions of some problem). They are both finite and I can easily work with them, for example plot them using matplotlib. However when I try to work with function f_3(k, \epsilon) which is a product of f_1 and f_2, i.e.

f_3 := f_1 f_2,

then this doesn’t work. For example if I want to plot f_3 I get the following error:

'''
snippet from my code that demonstrates the problem
'''
c = plot(fun1, title = 'fun 1')                    # No problem plotting fun1
plt.colorbar(c)
plt.show()

c = plot(fun2, title = 'fun 2')                    # No problem plotting fun2
plt.colorbar(c)
plt.show()

c = plot(fun1 * fun2, title = 'fun 1 * fun 2')     # Cannot plot fun1 * fun2. WHY??
plt.colorbar(c)
plt.show()

Traceback (most recent call last):
  File "/mnt/c/Users/juraj/Desktop/Master Thesis/code/finished-models/min_ex.py", line 142, in <module>
    c = plot(fun1 * fun2, title = 'fun 1 * fun 2')     # Cannot plot fun1 * fun2. WHY??
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py", line 436, in plot
    return _plot_matplotlib(object, mesh, kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py", line 302, in _plot_matplotlib
    return mplot_function(ax, obj, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py", line 152, in mplot_function
    return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)
  File "/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py", line 307, in tricontourf
    return TriContourSet(ax, *args, **kwargs)
  File "/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py", line 29, in __init__
    super().__init__(ax, *args, **kwargs)
  File "/usr/lib/python3/dist-packages/matplotlib/contour.py", line 812, in __init__
    kwargs = self._process_args(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py", line 45, in _process_args
    tri, z = self._contour_args(args, kwargs)
  File "/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py", line 72, in _contour_args
    raise ValueError('z array must not contain non-finite values '
ValueError: z array must not contain non-finite values within the triangulation

I can understand what the error says, i.e. that some value in f_3 is non-finite (at least that is how I understand it) but I just don’t get why this is since both f_1 and f_2 are finite (which I deduced by the fact that matplotlib had no issue plotting them).

I would very much appreciate if anyone would be able to look at this and try to help me find where the problem is. I am sorry for the lengthy minimal working example but I could not get the error when I tried with more simple code, the part where error occurs is at the very end of the code. This is the best I could do:

from fenics import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(4, 4)

# Create boundaries subdomains
walls   = 'x[1] < DOLFIN_EPS or x[1] > 1 - DOLFIN_EPS'
inflow  = 'x[0] < DOLFIN_EPS'
outflow = 'x[0] > 1 - DOLFIN_EPS' 

# Choose values for initial and boundary conditions 
u_init = 1 
k_init = 0.1   
e_init = 0.1 

# Create function spaces for all variables
V = VectorFunctionSpace(mesh, "CG", 2)       
Q = FunctionSpace(mesh, "CG", 1)                                        
K = FunctionSpace(mesh, "CG", 1)

# Create boundary conditions for u,k,e
bcu = [DirichletBC(V, Constant((0.0, 0.0)), walls)]
bck = [DirichletBC(K, Constant(0.0), walls)]
bce = []

# Create boundary conditions for pressure p
bcp_inflow = DirichletBC(Q, Constant(2.0), inflow)
bcp_outflow = DirichletBC(Q, Constant(0.0), outflow)
bcp = [bcp_inflow, bcp_outflow]

# Define trial functions
u = TrialFunction(V)
p = TrialFunction(Q)
k = TrialFunction(K)
e = TrialFunction(K)

# Define test functions
v = TestFunction(V)
q = TestFunction(Q)
phi = TestFunction(K)
psi = TestFunction(K)

# Define functions for quantities at previous time
u0 = interpolate(Constant((u_init, 0)), V)
k0 = interpolate(Constant(k_init), K)
e0 = interpolate(Constant(e_init), K)

# Define functions where solutions will be stored
u1 = Function(V)
p1 = Function(Q)
k1 = Function(K)
e1 = Function(K)

# Define constantss used in variational forms
dt = Constant(0.001)
nu = Constant(0.00181818)
 
# Calculate distance to the wall
dfun = Expression('0.5 - abs(0.5 - x[1])', degree = 1)
                  
# Define Min, Max function
def Max(a, b): 
    return (a+b+abs(a-b))/Constant(2)

# Define helping functions for k-eps model
Re_t = (1. / nu) * (k0**2 / e0)
Re_k = (1. / nu) * (sqrt(abs(k0)) * dfun)

f_nu = (1 - exp(- 0.0165 * Re_k))**2 * (1 + 20.5 / Re_t)
f_1  = 1 + (0.05 / f_nu)**3
f_2  = 1 - exp(- Re_t**2)
S_sq = 2 * inner(sym(nabla_grad(u1)), sym(nabla_grad(u1)))

# Define terms that appear in variational formulation for k-eps model
nut     = Max(0.09 * f_nu * (k0**2 / e0),    Constant(0.0))
prod_k  = Max(nut * S_sq,                    Constant(0.0))
react_k = Max(e0 / k0,                       Constant(0.0))
prod_e  = Max(1.44 * react_k * f_1 * prod_k, Constant(0.0))
react_e = Max(1.92 * react_k * f_2,          Constant(0.0))

# Create variational formulation for N-S (Chorin's splitting method)
F1  = (1. / dt) * dot(u - u0, v)*dx \
    + dot(dot(u0, nabla_grad(u)), v)*dx \
    + (nu + nut) * inner(grad(u), grad(v))*dx 

a_1 = lhs(F1); l_1 = rhs(F1)

a_2 = dot(grad(p), grad(q))*dx
l_2 = (-1. / dt) * dot(div(u1), q)*dx

a_3 = dot(u, v)*dx
l_3 = dot(u1, v)*dx - dt * dot(grad(p1), v)*dx

# Create variational formulation for k-eps
FK  = (1. / dt) * dot(k - k0, phi)*dx \
    + dot(dot(u1, nabla_grad(k)), phi)*dx \
    + (nu + nut / 1.0) * inner(grad(k), grad(phi))*dx \
    - dot(prod_k, phi)*dx + dot(react_k * k, phi)*dx 

FE  = (1. / dt) * dot(e - e0, psi)*dx \
    + dot(dot(u1, nabla_grad(e)), psi)*dx \
    + (nu + nut / 1.3) * inner(grad(e), grad(psi))*dx \
    - dot(prod_e, psi)*dx + dot(react_e * e, psi)*dx 

a_k = lhs(FK); l_k = rhs(FK)
a_e = lhs(FE); l_e = rhs(FE)

# main loop
iter_max = 3
tol = 1e-7

for iter in range(iter_max):

    solve(a_1 == l_1, u1, bcu)
    solve(a_2 == l_2, p1, bcp)
    solve(a_3 == l_3, u1, bcu)
    solve(a_k == l_k, k1, bck)
    solve(a_e == l_e, e1, bce)
    
    # Update all functions
    u0.assign(u1)
    k0.assign(k1)
    e0.assign(e1)

###################### PROBLEM HAPPENS HERE ######################
    
fun1 = (1. / f_nu)**2
fun2 = (1. / f_nu)

c = plot(fun1, title = 'fun 1')                    # No problem plotting fun1
plt.colorbar(c)
plt.show()

c = plot(fun2, title = 'fun 2')                    # No problem plotting fun2
plt.colorbar(c)
plt.show()

c = plot(fun1 * fun2, title = 'fun 1 * fun 2')     # Cannot plot fun1 * fun2. WHY??
plt.colorbar(c)
plt.show()

##################################################################

Thank you for reading and once again, I will appreciate any help I can get with this :smiley:.

2

There are a few things that happen under the hood here.

  1. Any plot of something that is not a dolfin.Function will be projected into P1 (of appropriate shape).
    Thus each of these plots is equivalent to:
fun1_p = project(fun1, V.sub(0).collapse())
fun2_p = project(fun2, V.sub(0).collapse())
sum_p = project(fun1*fun2, V.sub(0).collapse(), "mumps")

By inspecting each of these, you will observe that the values of each of these are way beyond machine precision (1e-32 to 1e-39) for fun1_p and (1e-17 to 1e-24) for fun2_p.

For a numerical point of view these values are 0, and will not be represented nicely in any assembly.

3

Hi @dokken, thank you for the fast reply :slight_smile:.
If I might, I have a few follow up questions.

Regarding what you wrote here, is it also true that if, for example, I use fun1 (the original fun1 that is not a dolfin.Function) inside a variational form (as a source term or diffusion coefficient) that it is also projected onto P1? If this is not true, then what happens under the hood when I use such UFL function inside variational form. I am asking this because this is exactly what I do in the minimal working example above, look for example at nut or prod_k in the minimal working example and how it is used in F1 or FK.

Secondly, are you aware if there is a better way of dealing with this? Is it ok to use such UFL functions inside variational formulation or should I always project onto some space. Is there maybe another way I am not seeing?

Thirdly, I am a bit confused about what “mumps” means, when I use it inside my code I get the following error:

TypeError: expected a (list of) <class 'dolfin.cpp.fem.DirichletBC'> as 'bcs' argument

could you elaborate?

Thank you for pointing this out, this makes sense to me why it’s not working. However, in my original script (I didn’t post it since it’s too long) the two functions f_1 and f_2 have much more reasonable values. I managed to change the MWO a little bit and now values of fun1 and fun2 are both around 1e+1. Yet, the error persists. Here is the new MWO where this is demonstrated (takes little longer to run, around 1min on my laptop):

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# Create mesh
height = 2.0; width = 2.0; nx = 4; ny = 512
def meshing(lx,ly, Nx,Ny):
    m = UnitSquareMesh(Nx, Ny)
    x = m.coordinates()

    #Refine on top and bottom sides
    x[:, 1] = 0.5 * ( np.arctan(2 * pi * (x[:,1] - 0.5)) / np.arctan(pi) + 1)

    #Scale
    x[:,0] = x[:,0]*lx
    x[:,1] = x[:,1]*ly
    return m

mesh = meshing(width, height, nx, ny)

# Create boundaries subdomains
walls   = 'x[1] < DOLFIN_EPS or x[1] > 2 - DOLFIN_EPS'
inflow  = 'x[0] < DOLFIN_EPS'
outflow = 'x[0] > 2 - DOLFIN_EPS' 
    
class Periodic(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)

    def map(self, x, y):
        y[0] = x[0] - 2
        y[1] = x[1]

periodic = Periodic()

# Choose values for initial and boundary conditions 
u_init = 20
k_init = 1.5  
e_init = 0.3

# Create function spaces for all variables
V = VectorFunctionSpace(mesh, "CG", 2, constrained_domain = periodic)       
Q = FunctionSpace(mesh, "CG", 1)                                        
K = FunctionSpace(mesh, "CG", 1, constrained_domain = periodic)         

# Create boundary conditions for u,k,e
bcu = [DirichletBC(V, Constant((0.0, 0.0)), walls)]
bck = [DirichletBC(K, Constant(0.0), walls)]
bce = []

# Create boundary conditions for pressure p
bcp_inflow = DirichletBC(Q, Constant(2.0), inflow)
bcp_outflow = DirichletBC(Q, Constant(0.0), outflow)
bcp = [bcp_inflow, bcp_outflow]

# Define trial functions
u = TrialFunction(V)
p = TrialFunction(Q)
k = TrialFunction(K)
e = TrialFunction(K)

# Define test functions
v = TestFunction(V)
q = TestFunction(Q)
phi = TestFunction(K)
psi = TestFunction(K)

# Define functions for quantities at previous time
u0 = interpolate(Constant((u_init, 0)), V)
k0 = interpolate(Constant(k_init), K)
e0 = interpolate(Constant(e_init), K)

# Define functions where solutions will be stored
u1 = Function(V)
p1 = Function(Q)
k1 = Function(K)
e1 = Function(K)

# Define constantss used in variational forms
dt = Constant(0.01)
nu = Constant(0.00181818)

# Calculate distance to the wall
dfun = Expression('1 - abs(1 - x[1])', degree = 1)

# Define Min, Max function
def Max(a, b): 
    return (a+b+abs(a-b))/Constant(2)

def Min(a, b): 
    return (a+b-abs(a-b))/Constant(2)

# Define helping functions for k-eps model
Re_t = (1. / nu) * (k0**2 / e0)
Re_k = (1. / nu) * (sqrt(abs(k0)) * dfun)

f_nu = Max(Min((1 - exp(- 0.0165 * Re_k))**2 * (1 + 20.5 / Re_t), Constant(1.0)), Constant(0.01116))
f_1  = 1 + (0.05 / f_nu)**3
f_2  = 1 - exp(- Re_t**2)
S_sq = 2 * inner(sym(nabla_grad(u1)), sym(nabla_grad(u1)))

# Define terms that appear in variational formulation for k-eps model
nut     = Max(0.09 * f_nu * (k0**2 / e0),    Constant(0.0))
prod_k  = Max(nut * S_sq,                    Constant(0.0))
react_k = Max(e0 / k0,                       Constant(0.0))
prod_e  = Max(1.44 * react_k * f_1 * prod_k, Constant(0.0))
react_e = Max(1.92 * react_k * f_2,          Constant(0.0))

# Create variational formulation for N-S (Chorin's splitting method)
F1  = (1. / dt) * dot(u - u0, v)*dx \
    + dot(dot(u0, nabla_grad(u)), v)*dx \
    + (nu + nut) * inner(grad(u), grad(v))*dx 

a_1 = lhs(F1); l_1 = rhs(F1)

a_2 = dot(grad(p), grad(q))*dx
l_2 = (-1. / dt) * dot(div(u1), q)*dx

a_3 = dot(u, v)*dx
l_3 = dot(u1, v)*dx - dt * dot(grad(p1), v)*dx

# Create variational formulation for k-eps
FK  = (1. / dt) * dot(k - k0, phi)*dx \
    + dot(dot(u1, nabla_grad(k)), phi)*dx \
    + (nu + nut / 1.0) * inner(grad(k), grad(phi))*dx \
    - dot(prod_k, phi)*dx + dot(react_k * k, phi)*dx 

FE  = (1. / dt) * dot(e - e0, psi)*dx \
    + dot(dot(u1, nabla_grad(e)), psi)*dx \
    + (nu + nut / 1.3) * inner(grad(e), grad(psi))*dx \
    - dot(prod_e, psi)*dx + dot(react_e * e, psi)*dx 

a_k = lhs(FK); l_k = rhs(FK)
a_e = lhs(FE); l_e = rhs(FE)

# main loop
iter_max = 88
tol = 1e-7

for iter in range(iter_max):

    solve(a_1 == l_1, u1, bcu)
    solve(a_2 == l_2, p1, bcp)
    solve(a_3 == l_3, u1, bcu)
    solve(a_k == l_k, k1, bck)
    solve(a_e == l_e, e1, bce)
    
    # Update all functions
    u0.assign(u1)
    k0.assign(k1)
    e0.assign(e1)

###################### PROBLEM HAPPENS HERE ######################
fun1 = f_nu
fun2 = k0**2 / e0

c = plot(fun1, title = 'fun 1')                    # No problem plotting fun1
plt.colorbar(c)
plt.show()

c = plot(fun2, title = 'fun 2')                    # No problem plotting fun2
plt.colorbar(c)
plt.show()

c = plot(fun1 * fun2, title = 'fun 1 * fun 2')     # Cannot plot fun1 * fun2. WHY??
plt.colorbar(c)
plt.show()
##################################################################

Lastly, I just want to clarify that the issue is not really about plotting the function f_3 = f_1 f_2. The issue is that function f_3 (which in reality is the function nut from the code) is used in the variational formulation as part of the diffusion coefficient in problems F1, FK and FE. However, after some time (88 iteration in this case) something breaks which causes f_3 to have non-finite values which is a huge problem since it appears in variational formulation.

PS: I just noticed a strange behaviour that even though nut from the code cannot be plotted (or used for anything), function prod_k = nut * S_sq can. This can also be seen by running the following without problem:

fun3 = S_sq
c = plot(fun1 * fun2 * fun3, title = 'fun 1 * fun 2 * fun3') 
plt.colorbar(c)                        # no problem plotting fun1 * fun2 * fun3
plt.show()                             # even though fun1 * fun2 cannot be plotted

Sorry for such a long reply, thank you again for responding :smiley:.

4

They are evaluated at quadrature points.

You do not need to project onto any space in the variational forms.

Just remove it. Typo on my end.

For the remainders of your question, there is a lot to digest here, and I do not have the bandwidth to do that at the correct instance. You are solving several PDEs, where nut depend on the previous solution e0. As you have a termination criteria of 88 iterations, do the plots look fine for any smaller subset of iterations?

I would strongly suggest that you add the three plots to your post as well, as one might be able to see a trend in them.

5

Thank you @dokken, your comments to my first 3 questions very very helpful :+1:.

Indeed nut depends on e0 as well as previous solution k0 (as do all the other quantities) an it is actually rather complicated function of those solutions. That’s why it’s been so difficult for me to find out where exactly the problem lies. As for your question, yes the plots look fine up until that point, it only breaks at the 88th iteration.

I am posting 5 plots here, first 3 are plots of f_1, f_2 and f_3 = f_1 f_2 from the last iteration where everything worked (i.e. iteration 87):
Fun1_working
Fun2_working
Fun1Fun2_working

The last 2 plots are of f_1, f_2 from iteration 88, i.e. the iteration where something goes wrong and f_3 = f_1 f_2 can no longer be plotted (or used effectively):
Fun1
Fun2

As you can see not much is different between the two iterations, that is why this seems really strange to me.

6

This looks weird to me. However, as your problem is rather complex, it is hard for me to give any further guidance. You would have to reduce the problem to a more tangible example for me have bandwidth to investigate further.