Incorrect Values Creating Expression From Other Expressions

Legacy DOLFIN
1

I am trying to create a more complicated expression from less complicated expressions following the examples I found here and elsewhere. I appear to be running into an issue with the keyword argument expression outputting the correct values functions of time. When I compare this “expression of expressions” to an explicit expression, they output different values. I have attached a MWE below demonstrating the problem below. I would like to know what I need to change so that they output the same value. Thank you for the help.

Edit: Upon further reflection, I think this is due to the interior expression not having it’s t updated when I update the t of the outer expression. Is there a way to pass the t-update through the outer expression to all interior expressions? Or do it have to update the interior expressions one by one, even if there are many of them, each time I change the t I want to evaluate the expression at?

FEniCS Version: 2019.2.0.13.dev0

from fenics import *
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, 'CG', 6)
a = Expression('cos(t)', t=0.0, degree=4)
test1 = Expression('cos(t)', t=0.0, degree=4)
test2 = Expression('a', a=a, t=0.0, degree=4)
test1.t = 0.1
test2.t = 0.1
t1 = interpolate(test1, V)
t2 = interpolate(test2, V)
print('t1(0.1) = ', t1(0.0))
print('t2(0.1) = ', t2(0.0))
2

You do have to update all of them yourself.

If the calculations are complicated, it may be easier to have them carried out by a symbolic math library, and then wrap just the final result in an expression.

from fenics import *
import sympy

t = sympy.var("t")
exp1 = sympy.sin(t)
exp2 = sympy.cos(t)
test1_sym = exp1 * exp2
test2_sym = exp1 + exp2
test1_cpp = sympy.printing.ccode(test1_sym)
test2_cpp = sympy.printing.ccode(test2_sym)


mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, 'CG', 6)
test1 = Expression(test1_cpp, t=0.0, degree=1)
test2 = Expression(test2_cpp, t=0.0, degree=1)
test1.t = 0.1
test2.t = 0.1
t1 = interpolate(test1, V)
t2 = interpolate(test2, V)
print('t1(0.1) = ', t1(0.0), "vs", test1_sym.subs({t: 0.1}))
print('t2(0.1) = ', t2(0.0), "vs", test2_sym.subs({t: 0.1}))

gives

t1(0.1) =  0.09933466539753065 vs 0.0993346653975306
t2(0.1) =  1.0948375819248541 vs 1.09483758192485
3

I will post a code later when I have the time, but the update is automatic when I use a Nonlinearproblem. After having solved for 2 unknows from a mixed element space, any expression I define involving those 2 variables are automatically up to date, there’s no need to invoke sympy or other libs.

4

That’s a different case. Here t is not a solution to another problem, it’s just a constant (the time, I guess)

5

Thank you! This helped tremendously, but now I am running into another issue using this approach. I can’t seem to write “a=a” inside the arguments for the Expression when using SymPy. This has never been a problem in FEniCS before. I also found out that if I change it to be “a=a0” inside the Expression, everything appears to work. I’m not sure what SymPy is doing in the background to break this functionality. It’s a minor issue, easily remedied, but not knowing what is going on concerns me and I would like to know the behavior so I avoid any issues in the future. I have attached an MWE demonstrating that a=a0 works, with f1, and a=a doesn’t, with f2.

from fenics import *
import sympy

t0 = 0.0
a = 1.0
a0 = 1.0
t_sym = sympy.var('t')
a_sym = sympy.var('a')
f_sym = -sympy.sin(t_sym) + a_sym
f_cpp = sympy.printing.ccode(f_sym)
print('f12_cpp = ', f_cpp)
f1 = Expression(f_cpp, t=t0, a=a0, degree=4)
f2 = Expression(f_cpp, t=t0, a=a, degree=4)
6

help(sympy.var) explains what is going on (arrows are my addition)

Help on function var in module sympy.core.symbol:

var(names, **args)
    Create symbols and (---->) inject them into the global namespace (<-----)

Simplest workaround is to swap the definitions of a_sym and a, so that a = 1.0 overwrites the variable created by sympy. Proper solution would be to look into sympy documentation on how to create a symbol without having it injected in the global namespace.

7

Thank you so much for the assistance! That all worked.