glwagner / dedales Goto Github PK

I want to write

closure = ConstantSmagorinsky(delta=0.01*Lx, C=0.13)
model = BoussinesqModel(closure=closure)

Can we do this @kburns?

Though our Chebyshev grid spacing is very fine near boundaries, we still probably want to implement wall models to compute fluxes for Dirichlet-type boundary conditions at some point.

@wenegrat and I talked about implementing a `split viscosity' method for stabilizing time-stepping.

Essentially, we would like to allow the user to specify some constant 'spliting viscosity' (or diffusivity) in the construction of an EddyViscosityClosure that is added to the left side of an equation and subtracted from the right side.

For example, introducing diffusivity splitting into the buoyancy equation produces

b_t + κ_split Δb + ... = ... + ∇ · ( κ_sgs - κ_split ∇b)

@wenegrat has found this technique allows larger time-steps and increases numerical stability.

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*ν.

@SandreOuza has found that dedaLES will run stably only if the ratio between the grid scale Δ and the Kolmogorov scale, defined as

ℓ = ∜(ν³ / ε) ,

where ν is viscosity and ε is turbulent dissipation, is below a certain number (perhaps around Δ/ℓ = 5). We should document this empirical rule. (Related to #19).

We might consider implementing spectral filters, or a spectral vanishing viscosity method to remove energy from small-scales and (potentially) stabilize our solutions. This is a form of LES which has received some attention in the literature:

https://www.cscamm.umd.edu/people/faculty/tadmor/spectral_viscosity/Kar2_LES_SV_JCP163_2000.pdf

Please send more references if they are obtained.

Can we implement functions / methods that operate on equations in spectral space in dedalus?

Would the use of compound grids mitigate the numerical stability problems we are experiencing?

How would we implement a "building-block" strategy for models (eg, BoussinesqModel = NavierStokes + Buoyancy Conservation, somehow)?

If the problem solves for a flow perturbation around a background field --- for example, if we solve for u, where the total flow field is U = u + u_background --- we may want the closure to act on the 'full' velocity field U rather than just u.

There are two ways to do this perhaps. One is to write the equations in terms of a non-u variable (u_prime, for example), and then to introduce a substitution:

problem.substitutions['u'] = "u_prime + u_background"

Another way is to allow the closure to accept keyword arguments; for example, changing the ConstantSmagorinsky definition to

class ConstantSmagorinsky(EddyViscosityClosure):
    __init__(..., u_background = "0", ...):
        self.u = f"u + {u_background})

The strain substitution, for example, would then be

def substitute_strainratetensor(self, problem):
        ...
        problem.substitutions['Sxx'] = f"{self.u}x"

Or something like that...

What is the best way to implement this functionality?

We need test cases to

1. Verify that the turbulence closures are implemented correctly (the math is right)
2. Verify that we can reproduce results from the literature
3. Test the accuracy of closures (and potential modifications thereof) relative to DNS

I am not sure how to do 1. with dedaLES. Perhaps we can write a test that, for example, computes Fx_sgs given some sample velocity field. I'm not exactly sure how to do this in dedalus. One common method for CFD codes is to compute a forcing so that

u_t + ... = ... + Fx_sgs + forcing = 0

For some chosen u(x, y, z). Then verifying that u does not change when the system is time-integrated provides a test that Fx_sgs is correct. This seems possible in dedaLES if we add forcing to our equation sets.

For 2 perhaps Stratified Couette flow a la Vreugdenhil and Taylor 2018 is a good place to start? Perhaps we can also find a non-stratified test case (isotropic homogeneous turbulence?) for constant Smagorinsky.

For non-literature cases, 3 is more of a research subject, but may be important since we plan to generalize closures somewhat (for example, to implement perturbations from background flow).

We will probably want to implement the "Deardorff" subgrid scale model, which uses turbulent kinetic energy as a prognostic variable to estimate eddy diffusivities.

It seems that the Deardorff model is often used for much coarser resolutions than other LES models, where 'coarseness' refers to the ratio between the grid scale and the Kolmogorov scale. For example, in the 'free convection' test case in Van Roekel et al 2018, the Kolmogorov scale (estimated by equating buoyancy flux with turbulent dissipation rate) is about 0.002 meters, while the grid scale is 0.5 meters --- a factor of 250 larger.

Anisotropic Minimum Dissipation, on the other hand, seems to be used for cases in which the grid scale is approximately 5-6 times the DNS grid scale (which is, hopefully, a rough guide for the Kolmogorov scale). This is true in Rozema et al 2015 (with a ratio of 64/384 in their section C) and Vreugdenhil and Taylor 2018 (with ratios 64/256, 128/512, and 128/768).

It would be interesting to have the Deardorff model implemented to understand the apparent discrepency between a 'Deardorff LES' (with grid to Kolmogorov ratios of 200 or larger) and a 'eddy viscosity LES' (with grid to Kolmogorov ratios of 4-6).

Some references for the Deardorff-style closure are:

• Pressel et al 2015
• Moeng 1984
• Deardorff 1980
• Sullivan et al 1994

The simplest presentation is in Pressel et al 2015. It's not clear whether the model described there is identical to the model implemented in Sullivan or Moeng. Perhaps a good place to start is to draw up some documentation that describes the model we would like to implement.

cc @SandreOuza

AMD closure currently returning NaN results

I started a branch for adding a quick setup.py script to make dedaLES at least locally installable using setup.py and/or pip. Along with this we probably want to pick a license for the code. I think this should probably be GPL3 given the reliance on Dedalus, which is GPL3.

Hi, I got an error on my centos machine when installing the dedaLED and testing it by python3 -m dedalus test.

The error is as follow:

E TypeError: 'numpy.float128' object cannot be interpreted as an integer

/opt/anaconda3/envs/dedalus/lib/python3.8/site-packages/scipy/special/_basic.py:2382: TypeError

I don't know why but I really install the dedaLES along the guide.