 # UFL expression with noninteger powers

Hi,

I have the following form for finding the motion of a dilute particle cloud, subject to an applied body force and air drag. The air drag expression requires one to take an non-integer power of u, i.e. power(Re_p(u),0.687)

``````def C_DtimesU(Re_p):
return  24.0/k_Re*(1.0 + 0.15*pow(Re_p,0.637))

def F_drag(u,c):
Re_p = k_Re*abs(u)
return -0.5*A_p*rho_f*C_DtimesU(Re_p)*c*u

F = (c*(u-u_1) + u*(c-c_1))*k_t*v_u*dx + (c*u*u).dx(0)*v_u*dx - \
F_drag(u,c)*k_mp*v_u*dx - c*F_applied*k_mp*v_u*dx + \
(c - c_1)*k_t*v_c*dx + (c*u).dx(0)*v_c*dx
``````

Is there a way of doing this within ufl? Could I make a c++ expression for the power function?

Alex

Do you mean taking non-integer powers in `Expression`s defined via C++ strings? You can use the C math library `pow` function in those, e.g.,

``````from dolfin import *
ex = Expression("pow(x,2)",degree=2)
print(assemble(ex*dx(domain=UnitIntervalMesh(1))))
``````

Hi yes something like that but it should take a function of u as an input i.e Re_p

I got it working using an equivalent expression using the exponential function:

``````ReLn = Constant(ln(0.687))
def C_DtimesU(Re_p):
return  24.0/k_Re*(1.0 + 0.15*exp(Re_p*ReLn))``````

Does the python ** operator not work?

``````Re_p**0.637
``````

Hi volkerk I think it was working in a previous version of fenics.

When I try solving the problem now I get nan values in the iteration statistics. I read somewhere in the new documentation that it might only work for integer power now.

According to the latest ufl manual (and my own scripts), the ** operator works, see p. 15, second line of the document:

Maybe you have to do

``````from ufl import *
``````

Hi thank you for the reply. I checked the UFL manual and know you can use the operator, however I’m not sure it works for non-integer powers. When I have the following imports in my file:

``````import matplotlib.pyplot as plt
from ufl import *
from fenics import *
import numpy as np
import time
``````

And use Re_p**0.687 I get the following in the iteration statistics:

`Newton iteration 46: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)`

any updates on this? I have the same problem

Hi gigi, I got it working using this equivalent expression:

``````ReLn = Constant(ln(0.687))
def C_DtimesU(Re_p):
return  24.0/k_Re*(1.0 + 0.15*exp(Re_p*ReLn))
``````

I had the same issue. However, I found out what was my mistake. Maybe you have the same mistake too. you can not have a negative number with non-integer power. Herein, I consider all of your parameter values are positive. So, Sometimes for numerical solvers maybe the solver estimates a negative value for your parameters which raises the issue. To solve this issue, you can easily use `abs` function for your parameters to fix the issue.

1 Like

Hi, I take the absolute value of the speed before it goes into the function as Re_p so that is taken care of. Maybe that part of UPL was fixed I wrote this work around a while ago