Unable to raise some Functions to non-integer powers

1

Hello,

I seem to be unable to raise exponential functions to non-integer values in FEniCS: they just return an array of nans.

I cannot seem to find any previous questions or documentation concerning this, so have come here.

For context, I have been attempting to implement a solution of non-linear diffusion-like equation, following https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/nonlinear-poisson/python/documentation.html.

I hope I haven’t made a silly mathematical error, but I am not sure what else to try.

from dolfin import *  # from fenics import * replaced by dolfin

# Setup:

DEG = 5
mesh = IntervalMesh(100, 0.0001, 20)
V = FunctionSpace(mesh, 'CG', DEG)

# First function, exponential, doesn't work
u_0 = Expression('exp(-pow(x[0], 2))', degree=DEG)
u_n = project(u_0, V)
powered=  pow(u_n, 2/3) 
print(project(powered, V).vector().get_local())  # A full array of NaNs

# Same function, exponential, integer power, works
u_0 = Expression('exp(-pow(x[0], 2))', degree=DEG)
u_n = project(u_0, V)
powered=  pow(u_n, 2) 
print(project(powered, V).vector().get_local())  # Working

# Linear, non-integer power, works
u_0 = Expression('x[0]', degree=DEG)
u_n = project(u_0, V)
powered=  pow(u_n, 2/3) 
print(project(powered, V).vector().get_local())  # Working

Many thanks.

$ dolfin-version
2019.1.0

Python running in a Jupyter notebook on Ubuntu 18.04.4

Python 3.6.9
[GCC 8.4.0] on linux

3

Looks like you’re failing to resolve a boundary layer with a coarse mesh in the projection. Either use more elements, or specify the exponent as a rational fraction so FFC can implicitly use optimisations. Cf.

from dolfin import *

DEG = 5
mesh = IntervalMesh(400, 0.0001, 20)
V = FunctionSpace(mesh, 'CG', DEG)

u_0 = Expression('exp(-pow(x[0], 2))', degree=DEG)
u_n = project(u_0, V)
powered=  pow(u_n, 2/3) 
print(project(powered, V).vector().get_local())

DEG = 5
mesh = IntervalMesh(100, 0.0001, 20)
V = FunctionSpace(mesh, 'CG', DEG)

u_0 = Expression('exp(-pow(x[0], 2))', degree=DEG)
u_n = project(u_0, V)
powered=  u_n**2/3
print(project(powered, V).vector().get_local())
1 Like
4

Thank you very much, I see how u_n no longer returns nans when you add more elements.

Though that works in the isolated example I gave, it does not seem to when actually implemented. Would you be able to help again?

I have simplified my code down to an extension of the non-linear poisson example from the documentation. I use 10,000 elements in an attempt to avoid the problem you mentioned, which I hope would be far more than enough for the interval [0,2].

The following code works for POWER=1, POWER=1.1, but not POWER=0.9 or POWER=0.5, where POWER is the exponent of u in the diffusion coefficient q(u).

import matplotlib.pyplot as plt
from dolfin import *

POWER = 0.9

T = 2.0 # final time
num_steps = 50  # number of time steps
dt = 0.0001#T / num_steps  # time step size

# Create mesh and define function space
DEG=5
nx=10000
mesh = IntervalMesh(nx, 0, 2)
V = FunctionSpace(mesh, 'CG', 5)

# Define initial value
u_0 = Expression('exp(-a*pow(x[0]-1, 2))', degree=DEG, a=5)
u_n = interpolate(u_0, V)

Diffusion coefficient

def q(u):
    return 1/(1+u**POWER)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Constant(0)
F = u*v*dx + dt*dot(q(u)*grad(u), grad(v))*dx - (u_n + dt*f)*v*dx

#####################################################
# Usual nonlinear problem formulation
du = TrialFunction(V)
a= derivative(F,u,du)# L is the weak form
problem= NonlinearVariationalProblem(F,u, J=a)
solver= NonlinearVariationalSolver(problem)

# Solver parameters
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-6
prm['newton_solver']['relative_tolerance'] = 1E-6
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0
#####################################################

t = 0
for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solver.solve()
    u_n.assign(u)
    if n in [1, 3, 7, 15, 30, 50]:
        plot(u, label="{:.2f}".format(t))
print("Success")

If I look at the log, I see that the residual immediately goes to nan for POWER=0.5

 Solving nonlinear variational problem.
    Newton iteration 0: r (abs) = 4.890e-03 (tol = 1.000e-06) r (rel) = 1.000e+00 (tol = 1.000e-06)
    Newton iteration 1: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
    ... [99 iterations removed]

Additionally, it also fails in the same way for q(u) = 1/(1+u*(u**0.5)), which could be useful information.

Many thanks again.

5

Your initial guess looks bad. You could probably use fewer elements, e.g. 1000 and just add the line:

u = Function(V)
u.interpolate(u_n)  # Set initial guess
v = TestFunction(V)
f = Constant(0)
F = u*v*dx + dt*dot(q(u)*grad(u), grad(v))*dx - (u_n + dt*f)*v*dx

to provide a suitable initial guess.

1 Like
6

I exactly have the same problem, do you have a solution for that?