Square Root and Natural Logarithm not Working

variational formulation
1

I’m running a FEniCS problem that involves solution of a diffusion-reaction pde across a meshed geometry. There are two pde’s being solved simultaneously with the independent variables ug and us, which represent concentrations in the gas and solid phases. We have been running the same code for a number of different reaction rate expressions. Most recently we are working with a rate expression that is described by the following function:

def rxn4(ug,us):
x = us/usmax
return (ug/c0) *
(A1 + B1exp(F1**2x) * (1.0 - erf(F1sqrt(x))) * (1/(exp((x-xm1)/s1)+1)) +
(A2 + B2exp(F2**2x) * (1.0 - erf(F2sqrt(x))) * (1 - 1/(exp((x-xm1)/s1)+1)) *
(1/(exp((x-xm2)/s2) + 1)) +
D*(1 - 1/(exp((x-xm2)/s2) + 1))

The independent variables are us and ug. A1, B1, A2, B2, F1, F2, D, s1, s2, xm1 and xm2 are constants. While attempting to use this function in our code, the exp() and erf() functions are working fine. However when I run with the sqrt() function I receive the following output:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.232e+02 (tol = 1.000e-09) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-09) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 2: r (abs) = nan (tol = 1.000e-09) r (rel) = nan (tol = 1.000e-09)
Newton iteration 3: r (abs) = -nan (tol = 1.000e-09) r (rel) = -nan (tol = 1.000e-09)

When I try to use the exp(ln(x)/2) I receive the same error.

I have tried to use sqrt(abs(x)) and exp(ln(abs(x))/2) to ensure that negative values are not passed to the sqrt and ln terms, but I receive the same output is received.

If I use exp(lnest(x)/2) to represent sqrt(x), where lnest(x) is a hyperbolic tangent summation series (Abramowitz & Stegun, eds. 1972 p. 68) function to estimate natural logarithm, this works fine. However at low values of x a large number of terms in the summation series are required to accurately represent ln(x) and this significantly increases run time. I would prefer to use either the sqrt or ln functions if possible.

Is there some reason that the exp() and erf() functions are working fine in this application while sqrt and ln functions are not working? Am I doing something wrong?

2

You might try adding a tiny epsilon to the argument, e.g.,

def safe_sqrt(x):
    return sqrt(x + DOLFIN_EPS)

I’ve run into problems before with x exactly equal to zero, where sqrt(x) is fine, but its symbolic derivative (as computed automatically when using the derivative function) involves a divide-by-zero.

2 Likes
3

Thank for this! Works fine.