again, "minor" trouble in detect.F90 & Re: [FLASH-BUGS] "minor" glitch in flatten.F90 (part 2)

[Tomasz Plewa]

Markus -

> Similar trouble as before, but now in detect.F90. Answer to the points you 
> raised follow below.
> 
> The problem in detect.F90 is:
> 
> line 129:
>     do i = 3, nzn5
>        if (scrch2(i-1)*scrch2(i+1) >= 0.e0) scrch3(i) = 0.e0
>     end do
> 
> in my case this flips the sign scrch2(i-1)*scrch2(i+1)=1.d-24 due to rounding 
> and produces quite significant differences in the solution, e.g. +-100 Pa 
> where the pressure in the problem ranges from 80000 Pa - 100000 Pa.

In essence this is pure arithmetics without any knowledge of
physics. If you question a 100 Pa difference produced by numerics then
perhaps your real question might be why 1.d-24 arises in your solution
and gains power to affect results. I do not think that in all cases
adding some tolerance is the cure, this one made me questioned how the
physical problem has to look like for such situation to take place.

> Additionally I noted:
> line 113:
>        if (a(i+1) - a(i-1) == 0.e0) then
>           temporary = small * smalla
>        else
> 
> I highly doubt that this will be true in a lot of cases. Is there any 
> rational for such a criteria, i.e. is it only here to avoid a div_by_zero ?

If the above condition is ever true then we are looking at density
field which is likely to be non-monotonic and no contact steepening
should take place. Moreover, monotonicity should completely flaten
density profile. So my guess is that we are only trying to avoid
division by zero from occuring.

> > I am generally against using floor-values in the code as they are (1)
> > problem-dependent, and (2) mask some real problems.
> 
> (2) I perfectly agree, this is exactly why I am following these up, wouldn't 
> waste days of my time if I would think, well it's just rounding error. From 
> my experience only 10% of "rounding errors" are rounding errors.
>
> 
> > 1.e-10 or 1.e+10 does not really mean much when taken out
> > of context. Hardwiring -10/+10 and -2/+2, as you suggested, may not be
> > the best solution. Tying those constant to epsiln may work better.
> 
> adding these though gives the uses some appreciation of specific trouble 
> makers and an easy way of ruling these out. In my case I would rather live 
> with a hardwired constant than a rounding error of 100.d-0 in the pressure... 
> epsilon unfortunately does not cover all rounding errors either, especially 
> if you have to set convergence criteria for iterative solutions where it may 
> not be necessary and or desirable to converge to machine precision. 
> Nevertheless, it is an excellent default for all hardwired constants.

I agree with your other. At the same time, I would also suggest to be
extra careful when we arrive to the point where numerics becomes a
dominant factor in our problem - perhaps we reached a limit when the
algorithm itself becomes sensitive to small perturbations.  Trying to
fix it is one thing, and great if possible (like in some of the cases
you mentioned), but it seems to me there exists a limit to that
approach as well.

Tomek
--
Tue, 13:30 CDT (18:30 GMT), Jul-22-2003
_______________________________________________________________________________

   Tomasz Plewa                                      www:   flash.uchicago.edu
   Computational Physics and Validation Group        email: tomek at uchicago.edu
   The ASCI Flash Center, The University of Chicago  phone: 773.834.3227
   5640 South Ellis, RI 475, Chicago, IL 60637       fax:   773.834.3230
_______________________________________________________________________________