Fortran round-off error with floating point operations

I have simple code, which flags nodes with in region enclosed by cylinder. On implementing the code, the result is mild tilt of the cylinder observed case with $\theta=90^{\circ}$.

The algorithm for checking any point inside arbitrarily oriented cylinder is as follows. Let $\vec{r}$ be the vector joining center $\vec{c}$ and arbitrary point $\vec{o}$, and $$\vec{r}=\vec{c}-\vec{o}.$$ For orientation vector $\vec{o}$, the projection of $\vec{r}$ on it is, $$u = \vec{r}\cdot\vec{o}.$$ Therefore, the perpendicular vector is, $$\vec{p}=\vec{r}-u\cdot\vec{o}.$$ For a cylinder of length $2l$ and radius $a$, check $$\vec{p}.\vec{p}<a^2 \hspace{1cm} \text{for } -l\leq u \leq l.$$

The actual issue: The above algorithm is implemented in Fortran. The code checks for points in Cartesian grid if inside the cylinder. Following being the test case: The cylinder makes an angle $\theta=90^{\circ}$ in the yz-plane with respect to y-axis. Therefore, the orientation vector $\vec{o}$ is (0, 1, 0).

Case 1: Orientation vector is assigned directly with $\vec{o}=(0.0,1.0,0.0)$. This results in perfect cylinder with $\theta=90^{\circ}.$

Case 2: Orientation vector is specified with intrinsic Fortran functions with double precision accuracy dsin and dcos with $\vec{o}=(0.0, \sin(\pi/2.0), \cos(\pi/2.0))$ with $\pi$ value assigned with more than 20 significant decimal points. The resulting cylinder results in mild tilt.

The highlighted region indicates the extra material due to tilt of cylinder with respect to Cartesian axes. I also tried architecture specific maximum precision "pi" value. This also results in same problem.

This shows like the actual angle made by cylinder is not $90^\circ$. Can anyone suggest valid solution for this problem. I need to use the inbuilt trigonometric functions for arbitrary angles and looking for accurate cell flagging method.

Note: All operations are performed with double precision accuracy.

• Though not related to the question, could you tell me how you made that nice picture? (Now I'm looking for the way to make such a picture for my study...) – roygvib Jun 18 '16 at 13:34
• Paraview and use color palette as "print" – Sathish Sanjeevi Jun 18 '16 at 13:41
• Could you be a little more precise about what you are drawing exactly? It would be much easier to help if there is a minimal reproducible example. Is $a$ an integer? Why does the cylinder look like a union of small cubes? What are the "nodes"? Are you doing something like testing whether points with integer coordinates belong to the cylinder? It could just be that your comparisons fall right on the boundary of the cylinder, which is ill-conditioned. – Kirill Jun 18 '16 at 23:30
• Those looking for the code for this could start at the cross-posted version -- stackoverflow.com/questions/37897677/fortran-round-off-errors. – High Performance Mark Jun 19 '16 at 9:58

It would be really nice to look at the code, compiler you are using and compiler options. However, I can suggest couple of things to test right away:

1. Instead of specifying 0.0 and 2.0 as in your text, use 0.0_dp and 2.0_dp, whenever you use any constants. That will make sure that double precision will be used for them. (dp=kind(0.0d0)) Moreover, I would say that multiplication by 0.5_dp is preferred to division by 2.0_dp. Say,

tmp_sin=dsin(pid*0.5_dp)
tmp_cos=dcos(pid*0.5_dp)
o=/0.0_dp,tmp_sin,tmp_cos/

2. Another thing to try, would be to create a double-precision variable pid_over_two that will have a value of pi/2 with >20 significant digits and try your code (with your compiler settings) on that.

3. Make sure that you do not compile this code with -no-prec-div option that allows for non-IEEE standard division.

4. Make sure that when you output the data for visualisation, you are using the correct format specifiers

These are the basic advices I can suggest without looking inside of the code. Also, it would be valuable to see, what output does your code (under your compiler/platform/settings) produce for:

dsin(pid/2.0)
dsin(pid/2.0_dp)
dsin(0.5_dp*pid)
dcos(pid/2.0)
dcos(pid/2.0_dp)
dcos(0.5_dp*pid)