[d3] Trying to understand taus practically

[c-white]

I've successfully run the Marti dynamo test, and everything appears to be working great. However, the actual problem I want to solve (still a dynamo on a ball) will have slightly different equations. I don't quite understand the part of the code that overwrites part of L with tau_columns, and so I don't know how I should modify it. Specifically:

1. When calling BC_rows(Nmax, ell, 9), is the 9 counting the first 9 equations (evolution + constraint, but somehow not the BCs proper that are enrolled later), or the 9 equations without condition="ntheta == 0", or something else?
2. How is it that N0-N2 are associated with u, N4-N6 with A, and N7 with T? What about N3 and N8?
3. Does the second index of L simply correspond to the fields in the order in which they are enrolled in the IVP?
4. Is the ell == 0 case always the same as ell != 0 but with the vector fields dropped from consideration?

[kburns]

Hi Chris. I wouldn't put too much effort into understanding the manually entered tau terms -- we have a new system that automates this programmatically and should be much easier to use and understand. Not all of the examples are updated to use this yet, but they should be soon. Here's one that's updated: test_new_marti_conv.py. You can see that you can directly call the LiftTau operator in the relevant equation, instead of having to track down the right entries.

[kburns]

Daniel is currently going through and updating most of the examples to use the new LiftTau structure (in the d3_examples branch). There are a few examples that still need some work, but these should be done and merged soon.

[kburns]

We now have narrative documentation on how to implement gauge conditions and tau terms directly in the problem formulation. Please feel free to reopen for further discussion if anything is unclear!

[hdrake]

General comments

The gauge conditions documentation is pretty clear, but I am still a bit confused about how to properly implement the tau method in the problem formulation.

It would be helpful to have more guidance as to how users should pick a basis to lift the tau terms to, and how that relates to the python methods arguments (e.g. in the 2D Poisson example). In both the documentation and the current d3 examples, the "lifted" basis is also sometimes modified by passing arguments a, b, k to the .clone_with method, without a clear explanation of why or what these parameter values mean (perhaps adding a Docstring to the method would help). I went through all of the examples scripts but it did not clarify much, except that the tau terms evidently need to be set differently for different problem types.

Specific questions

In this d3 notebook, I tried re-implementing @wenegrat's d2 notebook from his 2018 publication on bottom mixed layer instabilities (which is a similar problem to the Boccaletti et al. 2007 surface mixed layer instability example illustrated in the Eigentools documentation). I am able to replicate Wenegrat's LBVP results fine, but the EVP problem results are nonsensical (e.g. no positive imaginary eigenvalues when there should be). I've gone over the code many times to find the bug and it seems most likely the problem resides in how I implemented my tau terms.

## How should I be setting these?
lift_basis = zbasis.clone_with(a=1/2, b=1/2) # First derivative basis
lift = lambda A: d3.Lift(A, lift_basis, -1)

# First-order reductions
uz = dz(u) + lift(τ_uz)
vz = dz(v) + lift(τ_vz)
wz = dz(w) + lift(τ_wz)
bz = dz(b) + lift(τ_bz)

# setup problem
problem = d3.EVP(
    [u, v, w, b, p, # perturbation variables
     τ_u, τ_v, τ_w, τ_b, τ_p, τ_uz, τ_vz, τ_wz, τ_bz], # taus for enforcing boundary conditions
    eigenvalue=ω,
    namespace=locals()
)
problem.add_equation(# Cross-slope momentum
    'dt(u) + w*dz(U) + U*dx(u) + V*dy(u) - f*v*cosθ + dx(p)'
    '- b*sinθ - σ*dz(κ*uz)'
    '+ lift(τ_u) = 0'
)
problem.add_equation((# Along-slope momentum
    'dt(v) + w*dz(V) + U*dx(v) + V*dy(v) + f*(u*cosθ - w*sinθ) + dy(p)'
    '- σ*dz(κ*vz)'
    '+ lift(τ_v) = 0'
))
problem.add_equation((# Slope-normal momentum
    'dt(w) + U*dx(w) + V*dy(w) + f*v*sinθ + dz(p)'
    '- b*cosθ - σ*dz(κ*wz)'
    '+ lift(τ_w) = 0'
))
problem.add_equation((# Buoyancy
    'dt(b) + u*N2*sinθ + w*N2*cosθ + w*dz(B) + U*dx(b) + V*dy(b)'
    '- dz(κ*bz)'
    '+ lift(τ_b) = 0'
))
problem.add_equation('dx(u) + dy(v) + wz + τ_p = 0') # Continuity, modified with pressure tau

# Bottom BC
problem.add_equation('u(z=0) = 0')
problem.add_equation('v(z=0) = 0')
problem.add_equation('w(z=0) = 0')
problem.add_equation('dz(b)(z=0) = 0')

# Upper BC
problem.add_equation('dz(u)(z=H) = 0')
problem.add_equation('dz(v)(z=H) = 0')
problem.add_equation('w(z=H)*cosθ + u(z=H)*sinθ = 0') # vanishing vertical (not slope-normal) velocity
problem.add_equation('dz(b)(z=H) = 0')

problem.add_equation("integ(p) = 0") # Pressure gauge

solver = problem.build_solver()
solver.solve_dense(solver.subproblems[0])

Any thoughts or advice, @kburns?

[kburns]

Yeah EVPs can be more sensitive to the tau choices. But here I think there may be some other things. First, the variables in the EVP are implicitly multiplied by horizontal Fourier modes, so except when k=l=0, the variables have no "mean", so you actually don't need τ_p or to set the pressure gauge here.

However I'm still having trouble reproducing the Wenegrat 2018 results from that script even with Dedalus 2 -- it's possible that since the problem is formulated in meters and seconds, the matrix entry and NCC coefficients are so small that they're being zeroed by the hard cutoffs (a behavior which changed in 2020 in Dedalus 2). This is fixed by passing the ncc_cutoff=1e-10 and entry_cutoff=0 keywords to the problems in Dedalus v2 to recover the old (2018) defaults, but I'm still not quite seeing the right behavior...

[wenegrat]

I doubt this is the fundamental problem, but note @hdrake that the wavenumbers reported in the paper are in units of inverse meters, such that the fastest growing mode (in the code using rad/m) would be expected near l=1.3e-3.

That being said, the choice of l=2e-4 in your notebook now does look like it should fall above the low-wavenumber cutoff for instability.

Sorry this isn't of much help, I look forward to following what you find.

[hdrake]

the wavenumbers reported in the paper are in units of inverse meters Thanks, @wenegrat for clarifying—I thought that might be the case.

[hdrake]

Yeah EVPs can be more sensitive to the tau choices. But here I think there may be some other things. First, the variables in the EVP are implicitly multiplied by horizontal Fourier modes, so except when k=l=0, the variables have no "mean", so you actually don't need τ_p or to set the pressure gauge here.

Thanks—fixed that.

However I'm still having trouble reproducing the Wenegrat 2018 results from that script even with Dedalus 2 -- it's possible that since the problem is formulated in meters and seconds, the matrix entry and NCC coefficients are so small that they're being zeroed by the hard cutoffs (a behavior which changed in 2020 in Dedalus 2). This is fixed by passing the ncc_cutoff=1e-10 and entry_cutoff=0 keywords to the problems in Dedalus v2 to recover the old (2018) defaults, but I'm still not quite seeing the right behavior...

Setting ncc_cutoff=1e-10 and entry_cutoff=0 gives me (approximately) correct results. I'm not sure what issues you were having Keaton but your suggestions worked for me!

From my updated d3 notebook:

From Jacob's paper (using d2):

Thanks @kburns and @wenegrat!