Use strict inequality for search radius check?

[efaulhaber]

We currently only check if distance2 <= search_radius^2. However, this causes problems in some computations. For example, this here in the surface tension computation:

julia> support_radius = 0.02; distance = support_radius; (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25
ERROR: DomainError with -6.938893903907228e-18:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
     ⋮ internal @ Base.Math
 [2] ^(x::Float64, y::Float64)
   @ Base.Math ./math.jl:1206
Use `err` to retrieve the full stack trace.

This does not happen when distance < search_radius. The question is if there are any downsides to this strict inequality check. I don't see any.

Note that distance2 < search_radius^2 is not sufficient:

julia> 0.0004 * (1 - 0.7eps()) < 0.02^2
true

julia> sqrt(0.0004 * (1 - 0.7eps())) < 0.02
false

So we would have to do something like

distance2 < search_radius^2 * (1 - eps())

[LasNikas]

For example, this here in the surface tension computation:

Do you have a link for the code?

I mean, if this is an edge case, why not handle the edge case differently instead of changing the check to strict inequality?

[svchb]

Its currently handled differently.

    # The neighborhood search has an `<=` check, but for `distance == support_radius`
    # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`.
    # TODO Change this in the neighborhood search?
    # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19
    distance >= support_radius && return zero(pos_diff)

    distance <= 0.5 * support_radius && return zero(pos_diff)

    # Eq. 7
    A = 0.007 / support_radius^3.25 *
            (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25

    # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as
    # `m_b = rho_0 * volume` (Akinci boundary condition treatment)
    adhesion_force = -adhesion_coefficient * m_b * A * pos_diff / distance

    return adhesion_force

The question is only if the special handling is necessary when we can instead reduce the search area by eps().

[LasNikas]

Yes, I see. But:

Note that distance2 < search_radius^2 is not sufficient:

which means we loos a bit of performance when we calculate the sqrt(distance2) and then check distance < search_radius.

[svchb]

It should be noted that < was sufficient in all practical cases.

[efaulhaber]

which means we loos a bit of performance when we calculate the sqrt(distance2) and then check distance < search_radius.

As I said, the following should be sufficient:

distance2 < search_radius^2 * (1 - eps())

[LasNikas]

I can live with distance2 < search_radius^2 * (1 - eps()) with an explanatory comment. The weight and thus the impact of particles that are far away is negligible anyway...