UFL expression with noninteger powers

variational formulation
1

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

2

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))))
3

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))
4

Does the python ** operator not work?

Re_p**0.637
5

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.

6

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 *
7

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)

8

any updates on this? I have the same problem

9

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))
10

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
11

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