Using AMD makes solution go to NaN when solution is zero (or laminar?)

[sandreza]

AMD closure currently returning NaN results

[glwagner]

Can you provide a code snippet that reproduces this result?

[sandreza]

Here is a piece of code that reproduces the results:

The errors occur at iteration 2 where the Reynolds number is NaN (after the first time step). Note that changing to other closure schemes yield finite answers for the Reynolds number. Previously I had been getting an error that stated

/opt/dedalus/src/dedalus/dedalus/core/operators.py:995: RuntimeWarning: divide by zero encountered in power np.power(arg0.data, arg1.value, out.data)
/opt/dedalus/src/dedalus/dedalus/core/operators.py:758: RuntimeWarning: invalid value encountered in multiply np.multiply(arg0.data, arg1.data, out.data)

The error can be reproduced by running the script on a single core. Going to multiple cores the error is no longer printed.

Code below

---

import sys; sys.path.append("..")
import time, logging
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
from dedalus.extras import flow_tools
import dedaLES
logger = logging.getLogger(__name__)
# Fixed parameters
Ra = 2000
nx = ny = nz = 8 
Lx = Ly = 6.0               # Horizonal extent
Lz = 1.0                    # Vertical extent
Pr = 0.7                    # Prandtl number
Δb = 1.0                    # Buoyancy difference
# Calculated parameters
ν = np.sqrt(Δb*Pr*Lz**3/Ra) # Viscosity. ν = sqrt(Pr/Ra) with Lz=Δb=1
κ = ν/Pr                      # Thermal diffusivity 
# Construct model
closure = dedaLES.AnisotropicMinimumDissipation()
model = dedaLES.BoussinesqChannelFlow(Lx=Lx, Ly=Ly, Lz=Lz, nx=nx, ny=ny, nz=nz, ν=ν, κ=κ, Δb=Δb, closure=closure, nu=ν,V=Lx*Ly*Lz)
model.set_bc("no penetration", "top", "bottom")
model.set_bc("no slip", "top", "bottom")
model.problem.add_bc("right(b) = 0")
model.problem.add_bc("left(b) = Δb")
model.build_solver()
# Initial condition: unstable buoyancy grad + random perturbations
z = model.z
b0 = z/ Lz
model.set_fields(b=b0)
model.stop_at(iteration=3)
# Flow properties
flow = flow_tools.GlobalFlowProperty(model.solver)
flow.add_property("sqrt(u*u + v*v + w*w) / nu", name='Re_domain')
dt = 1e-4
for i in range(2):
    model.solver.step(dt)
logger.info('Max domain Re = %e' %flow.max("Re_domain"))

[glwagner]

Working on this issue on the branch AMDNaNs.

This file reproduces the error when the parameters a=0. The error disappears when a=10.

The error returns even when a=10 when b is not set initially.

[glwagner]

@SandreOuza and I have discussed how this problem is associated with the denominators of the subgrid viscosity expression and subgrid diffusivity expression going to zero.

@kburns any ideas how to solve this problem? I'm not sure how to implement any of the solutions I can think of (mainly involving if statements to catch the case that the denominator is 0) in dedalus.

@SandreOuza has suggested either setting a maximum value for the viscosities and diffusivities, or adding some small constant to the denominators of the viscosity and diffusivity. I dislike the idea of hard coding something with units, like a constant shift of tr_uij. Possibly, such a fudge factor can be a parameter that is initialized when the closure is initialized, with some sensible default. I'm just not exactly sure how we'd choose that sensible default. I don't know how to programmatically select a constant that is "small relative to tr_uij most of the time except when tr_uij=0".

Ideas?

[kburns]

It looks like for the sub grid shear viscosity and sub grid diffusivity, the resulting stress/flux should be zero when the denominators go to zero, so maybe the right way is to just short-circuit the computation of these terms and set them to zero at those points. At first glance, though, it's not clear that the buoyancy-induced sub grid stress also goes to zero when the strain rate goes to zero, so there may be some trouble there...

[glwagner]

@kburns I think short-circuiting makes most sense. How'd you implement that in dedalus?

I think the buoyancy-induced term might sort out too... naively, the buoyancy-induced subgrid stress is only zero for arbitrary buoyancy fields when wx=wy=wz=0. If u=v=0, then we also have tr_uij=0.

[kburns]

Just to document our discussion on this... one thought is that we might want to just leave this be, since it may only come up for initial conditions, and any fix based on short-circuiting would require work to check for zeros at every tilmestep. So maybe the best solution is to always start from slightly noisy initial conditions?

[glwagner]

Agreed.