Round-off error while using double precision variables

[Nagarada_Gadde__Srin]

I have a parallelised MPI code which solves non-linear PDEs (basically LES simulations), written in Fortran 90. In a particular section of the code I have to multiply double precision variables with an other double precision variable.

        u(iwrap_p(:),j,k) = alpha(:) * u(iwrap_p(:),j,k) + beta(:) * u(iwrap_p(:),j,k)

where alpha and beta are double precision variables whose values vary between 0. and 1.0, u is the velocity,and the sum of alpha and beta at any position is 1.0. The output of this operation should be equal to u(iwrap_p(:), j, k). But, performing this operation i.e. multiplication of the velocity with alpha and beta introduces round-off errors which keep on growing with time and finally the solution becomes non-physical.

Alpha and beta are calculated by the following, 

   beta = 0.5_rprec * ( 1.0_rprec - cos (pi * real (i - istart, rprec)  &
        / real(iend - istart, rprec)) )

    alpha = 1._rprec - beta 

where 'rprec' is the  user defined double precision datatype. 

I am using the following compiler arguments

  FFLAGS = -O3 -fp-model source

I have tried compiling with and without optimisations, but nothing seem to work. 

The error disappears when alpha and beta are declared as single precision variables. What might be the cause of this error? 

[jimdempseyatthecove]

Try:

-fp-model precise -fimf-precision=high

Also assure that pi is correct for double precision.

Jim Dempsey

[mecej4]

Sorry, something got lost in the narrative. If alpha(:) + beta(:) = 1, the calculation of u can be skipped completely, because it accomplishes nothing (other than to put in some round-off error).

[Nagarada_Gadde__Srin]

@mecej4 I am using the operation just as a sanity check. The actual code is the following: 

u(iwrap_p(:),j,k) = alpha(:) * u(iwrap_p(:),j,k) + beta(:) * u_previous_domain(iwrap_p(:),j,k)

where u_previous_domain is the velocity from the precursor domain which is being copied to the next computational domain. It is basically the addition of a fraction of velocity from the precursor domain to the consumer domain which utilises the velocity. 

@jimdempseyatthecove I will try it and let you know if it works out.  

[jimdempseyatthecove]

Your statement in #4 appears to be an Adams integration adjustment of an array (indexed by an index array) where u was previously advanced using Euler integration. If so, why didn't you write u once using Adams integration? Sketch:

velocity(:) = velocity(:) + (alpha(:) * acceleration(:) + beta(:) * prior_acceleration(:)) * delta_t

Typically alpha and beta do not vary.

Jim Dempsey

[Nagarada_Gadde__Srin]

@jimdempseyatthecove Your suggestion in #2 did not work. The error still keeps growing. 
Sorry about the lack of information in my question, perhaps I over-simplified it. I am not doing any integration. The simulation has two computational domains, one produces atmospheric boundary layer flow by solving Navier-Stokes equation (precursor domain), the second domain has a wind turbine (consumer domain); velocity from the first domain has to be given as input to the turbine domain. In the absence of the turbine, velocities in both the domain have to be identical, even after copying. The following operation does the copying:  

u(iwrap_p(:),j,k) = alpha(:) * u(iwrap_p(:),j,k) + beta(:) * u_precursor_domain(iwrap_p(:),j,k)

Fractions alpha and beta correspond to the coefficients which slowly set the velocity in the consumer domain to the velocity in the precursor domain. Alpha varies from 1 to 0 and beta varies from 0 to 1 making the velocity at the exit of the consumer domain equal to the velocity at the exit of the precursor domain. For this test case as there is no turbine, u(iwrap_p(:), j, k) and u_precursor_domain(iwrap_p(:), j, k) are identical. 
I might be expecting too much after using alpha and beta which obviously introduces round-off errors. But, I don't see the error when I use a single precision or quadruple precision for alpha, beta and the calculation of those coefficients. The error renders the velocity non-physical. I get perfect, physically meaningful solutions for both the precursor and consumer domains if I use single or quadruple precision for the variables. 

Here are the typical values of alpha and beta: 

Alpha                                           Beta 
0.933012701892219                               6.698729810778065E-002
0.750000000000000                               0.250000000000000 
0.500000000000000                               0.500000000000000
0.250000000000000                               0.750000000000000
6.698729810778059E-002                          0.933012701892219
0.000000000000000E+000                          1.00000000000000
0.000000000000000E+000                          1.00000000000000
0.000000000000000E+000                          1.00000000000000
0.000000000000000E+000                          1.00000000000000   
                   

Output from one of the operations where the coefficients are defined with double precision. 

 u : Before           22.3470033697708

 u:  After              22.3470033697707

[jimdempseyatthecove]

I haven't run the following in Fortran so I cannot tell from your typical values of alpha and beta are correct. Incorporating those numbers into Excell shows line 06 above as producing a sum that is 1 lsb less than 1.0. If this tends to be the case (error rounds down), then I would check the code that generates and then force the condition that alpha(n) + beta(n) == 1.0. This may be facilitated by manipulating alpha(n) and/or beta(n) in a manner that statically is neutral in balancing the factors. (IOW tweak using the +/- EPSILON(alpha(n)) and/or +/- EPSILON(beta(n)) to produce 1.0 exactly.