Mitigating near-zero eddy viscosities

[glwagner]

We need a softmax:

1. to ensure an accurate spectral approximation of the AMD viscosity, and
2. potentially to prevent any near-zero eddy viscosity (in either AMD or Constant Smagorinsky) from destabilizing the solution

This may help address some of the problems discussed in #20.

@kburns and I talked about this a while ago and we did not settle on a solution, though I cannot remember why.

I propose that we add parameters either to each closure individually or to EddyViscosityClosure, perhaps called ν_soft and κ_soft (better names welcome). We could then write (either in add_substitutions_subgrid_stress in EddyViscosityClosure or in add_substitutions in each closure individually):

if ν_soft is not None:
    problem.substitutions['softstep(x, a, b)'] = "0.5 * (1 + tanh((x-a)/b))"
    problem.substitutions['softmax(x, a, b)'] = "x * softstep(x, a, b)"

    problem.substitutions['ν_sgs_pred'] = # physics 
    problem.substitutions['ν_sgs'] = "softmax(ν_sgs_pred, ν_soft, ν_soft)"

(And similarly for κ_sgs.)

As written above, a in softmax is a displacement and b is a width. The appropriate relative values of a and b can be debated --- perhaps a more conservative choice is a=2*b? We can also add two parameters for the softmax to __init__; maybe with a sensible default. That'd give users full flexibility.

We should also change zero_max to hardmax in AMD?

A good value for ν_soft is probably the molecular ν, or perhaps 2*ν.

[kburns]

I think that particular softmax function can still go negative except for large values of a/b. An alternative is softmax(x, b) = b * np.log(1 + np.exp(x/b)), but this can have numerical issues for x/w > a few hundred.

[glwagner]

Hm!

We could make softmax a kwarg that the user may change if they desire?

w=b in what you've written, right?

An alternative way of writing your function is

softmax(x, a) = x + a*np.log(1 + np.exp(-x/a))

would that have fewer issues for large x/a?

[kburns]

Yeah sorry, w=b (fixed). Sure we could make it an option. The issue is the exponential can underflow or overflow if abs(x/a) is out of the exponent range of doubles.

[kburns]

Here we go: softmax(x, a) = a * np.logaddexp(0, x/a).

[glwagner]

Nice find!

If you think that's better than x + a*np.log(1 + np.exp(-x/a)), I'm good with it.

It seems that the parameter a should be a "small" fraction of the molecular viscosity or diffusivity: the error at x=0 is 0.7a, which is a sizeable fraction of a? So with a=nu/100, you'll incurr an effective error in your molecular diffusivity of 0.7% everywhere; perhaps acceptable. Or do we think that creates too-sharp a cut-off? The more turbulent the flow, the larger value of a you can accommodate.

[glwagner]

Closed by #23.