FEniCS Project

Controlling the precision

Asked by mwelland

Hi all,

I am having a problem that I think is related to precision but I can't resolve it by changing the form_compiler-precision option. An example of my problem is below.

I am setting a variable u that is always 0<= u <=1 . I define another variable p(u) which is also 0<= p <=1 for the same range of u values.
I test the values of u and p by projecting a function including the ln of the values. u passes but p fails *unless, you insert an abs in the test of p (as below. Remove it to see the NaN error I get).

It looks to me like a problem with precision / truncation error but I'm not sure how to fix it. Does anyone have any advice please?

Thanks,
Mike

~~~~~~~~~~~~~~~
from dolfin import *

parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["precision"] = 300

mesh = UnitSquare(50,2)
V = FunctionSpace(mesh, "Lagrange", 1) # Order 1 function space

u = Function(V)

uIC = Expression(".5*(1-tanh((x[0]-.5)/.05))"); # Initial function. Always between 0 and 1

u.interpolate(uIC) # Set values
p = pow(u,3)*(6*pow(u,2)-15*u+10);
up = project(p,V)

test1 = project(ln(1e15*(1-u)+u),V) # Will only work if 0 <= u <= 1
test2 = project(ln(1e15*abs(1-p)+p),V) # Will only work if 0 <= p <= 1
~~~~~~~~~~~~~~~~~~~~~~~

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Anders Logg
Solved:
Last query:
Last reply:
Revision history for this message
Best Anders Logg (logg) said : 		#1

On Tue, Nov 06, 2012 at 03:55:59PM -0000, mwelland wrote:
> New question #213461 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/213461
>
> Hi all,
>
> I am having a problem that I think is related to precision but I can't resolve it by changing the form_compiler-precision option. An example of my problem is below.
>
> I am setting a variable u that is always 0<= u <=1 . I define another variable p(u) which is also 0<= p <=1 for the same range of u values.
> I test the values of u and p by projecting a function including the ln of the values. u passes but p fails *unless, you insert an abs in the test of p (as below. Remove it to see the NaN error I get).
>
> It looks to me like a problem with precision / truncation error but I'm not sure how to fix it. Does anyone have any advice please?
>
> Thanks,
> Mike
>
> ~~~~~~~~~~~~~~~
> from dolfin import *
>
> parameters["form_compiler"]["representation"] = "quadrature"
> parameters["form_compiler"]["precision"] = 300

I don't think setting it higher than 16 makes any difference. :-)

> mesh = UnitSquare(50,2)
> V = FunctionSpace(mesh, "Lagrange", 1) # Order 1 function space
>
> u = Function(V)
>
> uIC = Expression(".5*(1-tanh((x[0]-.5)/.05))"); # Initial function. Always between 0 and 1
>
> u.interpolate(uIC) # Set values
> p = pow(u,3)*(6*pow(u,2)-15*u+10);
> up = project(p,V)
>
> test1 = project(ln(1e15*(1-u)+u),V) # Will only work if 0 <= u <= 1
> test2 = project(ln(1e15*abs(1-p)+p),V) # Will only work if 0 <= p <= 1

Why the factor 1e15? It seems bound to give you trouble.

The projection of your expression will lead to evaluation at
quadrature or nodal points. Some of those may lie close to an x-value
where function value is 0 and so it's likely round-off errors may
cause problems. My advise would be to add some small eps to your uIC
expression to keep it safely away from zero, or live with the abs.

--
Anders

Revision history for this message
mwelland (akfypznt56ld2kjywwchehm8-mike-nw2wga2s5w6wnovqwws3u9v1) said : 		#2

Ah, okay...

Yeah, the large factor is a pain but unfortunately necessary (unless there is a nice way to take it out of the ln...). I had some success defining
p1m = -(u-1.0)**3*(6.0*(u-1.0)**2+15.0*(u-1.0)+10.0)

and then
test2 = project(ln(1e15*p1m+p),V)

but not a really huge difference.

Revision history for this message
mwelland (akfypznt56ld2kjywwchehm8-mike-nw2wga2s5w6wnovqwws3u9v1) said : 		#3

Thanks Anders Logg, that solved my question.

Revision history for this message
Anders Logg (logg) said : 		#4

On Wed, Nov 07, 2012 at 07:41:16PM -0000, mwelland wrote:
> Question #213461 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/213461
>
> Status: Answered => Solved
>
> mwelland confirmed that the question is solved:
> Ah, okay...
>
> Yeah, the large factor is a pain but unfortunately necessary (unless there is a nice way to take it out of the ln...). I had some success defining
> p1m = -(u-1.0)**3*(6.0*(u-1.0)**2+15.0*(u-1.0)+10.0)
>
> and then
> test2 = project(ln(1e15*p1m+p),V)
>
> but not a really huge difference.

Let C = 1e15. Then your expression is

  ln(C*(1 - x) + x)
= ln(C*(1 - x + x/C))
~ ln(C*(1 - x + eps))
= ln(C) + ln(1 - x + eps)

where eps is some small number you choose to avoid negative values in
1 - x. The value of eps doesn't really matter as long as you are not
close to x = 1. And I assume you won't need to be as close as 1/1e15?

--
Anders

Revision history for this message
mwelland (akfypznt56ld2kjywwchehm8-mike-nw2wga2s5w6wnovqwws3u9v1) said : 		#5

I tried that in the example problem and it does fix it without the need of the abs(), so thank you for that. Later in my main model though I end up taking the derivative of that expression and removing u from the last term does impact it unfortunately.

Thank you for your help. In the end, I am trying to reformulate the p(u) expression so that the model is more robust against the numerical error.