Division by zero in zero function case

[Arpit-Babbar]

In the smoothness indicator function, we can encounter division by zero in computing

energy = max((total_energy - total_energy_clip1) / total_energy,
                 (total_energy_clip1 - total_energy_clip2) / total_energy_clip1)

This will never occur if the indicator variable comprises of positive quantity but indicator variable may not always be positive (example, testing a jump discontinuity in linear advection or Burger's equation)

One way to handle this is by putting an if in case total_energy is zero (currently in PR, will extend to other files if that is okay). The better way is

tolE = eps(total_energy)
energy = max((total_energy - total_energy_clip1) / (total_energy + eps),
                 (total_energy_clip1 - total_energy_clip2) / (total_energy_clip1 + eps))

but this will technically change all other cases. Another option which won't interfere existing code is

tolE = eps(total_energy)
energy = max((total_energy - total_energy_clip1) / max(total_energy, eps),
                 (total_energy_clip1 - total_energy_clip2) / max(total_energy_clip1, eps))

All of these options do fix the issue in my test case of interest (1-D multiple profile).

[ranocha]

Thanks! Please find a comment below. Could you please also add a test for this?

[sloede]

I am wondering though if this is the right strategy. The indicator was designed to work with an "energy-like" quantity, such as the total energy. If you are using a different quantity that is not strictly positive, I am not sure if it still makes sense to use the indicator.

Wouldn't it be better to fix the positivity in the variable selection function? E.g., if you have a function that extracts the jump, can you just make it abs(jump) to always get a non-negative value?

[ranocha]

I am wondering though if this is the right strategy. The indicator was designed to work with an "energy-like" quantity, such as the total energy. If you are using a different quantity that is not strictly positive, I am not sure if it still makes sense to use the indicator.

Wouldn't it be better to fix the positivity in the variable selection function? E.g., if you have a function that extracts the jump, can you just make it abs(jump) to always get a non-negative value?

I don't understand what you say here. I guess this is meant for something like an identically vanishing solution, where you just have vanishing energy.

[ranocha]

Thanks! We need to wait for a hotfix of some CI issues caused by upstream changes.

[sloede]

I don't understand what you say here. I guess this is meant for something like an identically vanishing solution, where you just have vanishing energy.

I think I misunderstood something from the description of the PR, so what I wrote before does not really make that much sense, sorry.

I am still wondering: Does it really make sense to test for equality to zero in the general case, and what is the performance impact of that change? I have rarely encountered a situation where testing for zero makes sense in practice for floating point numbers. AFAICT, this will only happen for very special initial conditions (as otherwise it will never hit zero exactly), thus I wonder if it would make more sense to modify the initial conditions instead?

[Arpit-Babbar]

I don't understand what you say here. I guess this is meant for something like an identically vanishing solution, where you just have vanishing energy.

I think I misunderstood something from the description of the PR, so what I wrote before does not really make that much sense, sorry.

I am still wondering: Does it really make sense to test for equality to zero in the general case, and what is the performance impact of that change? I have rarely encountered a situation where testing for zero makes sense in practice for floating point numbers. AFAICT, this will only happen for very special initial conditions (as otherwise it will never hit zero exactly), thus I wonder if it would make more sense to modify the initial conditions instead?

Modifying the initial conditions is also fine but the library should give the user the right error when they do use a vanishing initial condition.

The previous works of Persson-Peraire (2006) and Klockner Et Al (2011) test such indicators on initial conditions that vanish at some points. So, I might naively expect Trixi's indicator to also work out of the box.

[sloede]

Modifying the initial conditions is also fine but the library should give the user the right error when they do use a vanishing initial condition.

The previous works of Persson-Peraire (2006) and Klockner Et Al (2011) test such indicators on initial conditions that vanish at some points. So, I might naively expect Trixi's indicator to also work out of the box.

You are right, the current situation is not satisfying insofar it leaves the user with weird errors and no documentation on this issue. I am just trying to figure out what the best solution is for everyone in Trixi.jl. Maybe it is indeed the approach chosen in this PR 🙂

[sloede]

Overall, this LGTM. My only concern is that this negatively affects performance, but since there are already a number of if statements in that loop, and some other expensive operations, I guess it's OK 🙂