Does the ordering of the terms in the linear form matter? [Solution: increasing the precision of the forms]

Question

Good morning!

I have the following linear form defined:

V = FunctionSpace(mesh, "DG", 1)
W = MixedFunctionSpace([DG]*n)
...
L = ... # some terms
L += 1/sqrt(4*PI)*Sigma_a[0]*ac*T**4*(Cv/(Cv + dt*Sigma_a[0]*4*ac*T**3))*v[0]*dx

Where:
- ac and Cv are some constants,
- v[0] is my first test function
- T is a function in V
- Sigma is defined as:

class MyXS_a(Expression):
    def eval(self, value, x):
        value[0] = ...

    def value_shape(self):
        return (1,)

Sigma_a = MyXS_a()

The code thus written runs but the results are wrong.
However, if I simply change the place of the first Sigma_a[0] term, namely if I replace:

L += 1/sqrt(4*PI)*Sigma_a[0]*ac*T**4*(Cv/(Cv + dt*Sigma_a[0]*4*ac*T**3))*v[0]*dx

with:

L += 1/sqrt(4*PI)*ac*T**4*(Cv/(Cv + dt*Sigma_a[0]*4*ac*T**3))*Sigma_a[0]*v[0]*dx

it works... Does it make any sense? In my more complicated code, it is really annoying because depending on where I add this Sigma_a[0] term, I get completely different answers...

I looked at the manual but it does not seem to mention any rules for the ordering of the terms of the forms.

Does anyone know what the problem with the first expression is?

Thanks a lot!
Vincent

What happens if you wrap the divisor in parenthesis and make 1 a double?

(1.0/sqrt(4.0*PI))

vs

1/sqrt(4*PI) — Charles, Jul 1, 2014
Unfortunately, it is not the reason but thanks for the suggestion! — V_L, Jul 1, 2014

mwelland · Answer 1 · Jul 1, 2014

Best answer

Could it be a precision error?
Try adding the following to the beginning of your code:

parameters["form_compiler"]["precision"] = 100

and see if it makes a difference. Also, are you using optimization? Maybe turn that off and see if it makes a difference...

answered Jul 1, 2014 by mwelland FEniCS User (8,410 points)
selected Jul 1, 2014 by V_L

It fixed it! I am really impressed!

What did this 'magic' line do exactly? Increase the precision of the forms when they are compiled? Not sure to understand what it concretely means.

In any case, thank you so much! I would not have found that myself :)

commented Jul 1, 2014 by V_L FEniCS User (4,440 points)

As far as I understand it, that is exactly what it does. For a more concrete answer maybe the developers can chime in?
Glad it helped anyway! BTW: I was told 16 was the max that affects it:
[https://answers.launchpad.net/fenics/+question/213461][1]

commented Jul 2, 2014 by mwelland FEniCS User (8,410 points)

Incidentally, the fact that you are finding this sort of error might mean you need to scale your variables (or functions) to prevent this sort of thing from happening. If you dissect that expression term by term (following bedmas of course) I think you'd find the problem.

commented Jul 2, 2014 by mwelland FEniCS User (8,410 points)

-f precision=n
Set the number of significant digits to n in the generated code. The default
value of n is 15.

commented Jul 2, 2014 by Charles FEniCS User (4,220 points)

Does the ordering of the terms in the linear form matter? [Solution: increasing the precision of the forms]

1 Answer