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

Comments

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)

---

Unfortunately, it is not the reason but thanks for the suggestion!

Answer 1

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

Comments

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

---

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]

---

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.

---

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