 |
 |
Aleksandar Donev wrote: > > Hello, > > Following my question about what <0.0 means for real numbers, here is a > little more ifnormation and another plea for help from those that can help. > > Here is a segment of code in a simulation which calculates the collision > time of two expanding spheres moving along linear paths. It is simply a > quadratic equation that needs to be solved, and the answer is C / > (-B+Sqrt(B**2-A*C)) if the two spheres collide (i.e. if this result is > nonnegative and exists). I need to be very careful that this gives a > very accurate answer, and tricks like replacing ">0" with > ">EPSILON(1.0)" have consistently lead to the simulation breaking down > (i.e. a misprediction of some collision), to which the simulation is > very sensitive, since it is an "event driven" molecular dynamics which > relies on proper ordering and prediction of the sequence of sphere > collisions. > > The code given below, which is a translation of the mathematical breakup > of the cases that can occur (the case C<=0 means the spheres are already > touching or overlapping and is a by-product of finite accuracy), runs > fine if I enable full IEEE handling: > > A = SUM (v_ij**2) - expansion_rate ** 2 I'm no expert here, but that's never stopped me before. Here you are letting the compiler do the sums in whatever order and then wanting the final accuracy to be smaller than epsilon. Is that reasonable? Wouldn't you be better off doing these sums in double precision (or double whatever r_wp is) to try to control the roundoff? Or maybe sort the values before doing the sums? Or maybe trying Kahan's trick of doing the additions to minimize round-off problems? It just seems odd to me to expect essentially perfect results from this longish series of floating point calculations. Dick Hendrickson > B = DOT_PRODUCT (r_ij, v_ij) - expansion_rate * current_diameter > C = SUM (r_ij**2) - current_diameter ** 2 > IF (((B <= 0.0_r_wp) .OR. (A < 0.0_r_wp)) .AND. (B**2-A*C >= > 0.0_r_wp)) THEN > IF (C <= 0.0_r_wp) THEN > IF (C <-Max(100.0_r_wp*EPSILON(1.0_r_wp), 1E-5_r_wp)) THEN > WRITE (ERROR_UNIT,*) "ERROR:", "Overlapping spheres ", > particle, " and ", partner, " detected at time ", & > & current_time, " at ", lsd_packing%lsd_time, ". > Overlap=", (Sqrt(SUM(r_ij**2))-current_diameter) / & > & current_diameter > lsd_in_error = .TRUE. > lsd_problem = (/ particle, real_partner /) > END IF > collision_time = current_time > ELSE IF ((B < 0.0_r_wp) .OR. (A < 0.0_r_wp)) THEN > collision_time = current_time + C / (-B+Sqrt(B**2-A*C)) > ELSE > WRITE (OUTPUT_UNIT,*) "Suspicios collision of particles > ", particle, " and ", partner, " not expected to collide: & > &A, B, C=", A, B, C > collision_time = HUGE (1.0_r_wp) > END IF > ELSE > collision_time = HUGE (1.0_r_wp) > END IF > > However, as NAG is telling me, suspicious IEEE operations are happening: > > Warning: Floating divide by zero occurred > Warning: Floating underflow occurred > > Debugging shows that the problem is at: > > Program received signal SIGFPE, Arithmetic exception. > 661 collision_time = current_time + C / > (-B+Sqrt(B**2-A*C)) > > Unfortunately I could not get the debugger (dbx90) to give me the values > of A, B and C, and so I am not sure what is happenning here. Other > compilers never complain about anything even if I use their switches to > enable IEEE signal trapping. > > How should I write something like this in Fortran 90? Or is using the > IEEE TR module and its facilities the only way to go? > > Thanks, > Aleksandar > > -- > __________________________________ > Aleksandar Donev > Complex Materials Theory Group (http://cherrypit.princeton.edu/) > Princeton Materials Institute & Program in Applied and Computational Mathematics > @ Princeton University > Address: > 419 Bowen Hall, 70 Prospect Avenue > Princeton University > Princeton, NJ 08540-5211 > E-mail: [log in to unmask] > WWW: http://atom.princeton.edu/donev > Phone: (609) 258-2775 > Fax: (609) 258-6878 > __________________________________