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?

Thank you for your time,
Alex

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

from dolfin import *
ex = Expression("pow(x[0],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:
https://readthedocs.org/projects/fenics-ufl/downloads/pdf/latest/

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)