Fortran round-off error with floating point operations

Question

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.

Answer 1

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)