Problems with Chebyshev grid

[Kesava Prasad]

Dear Deadalus users,

I am testing a working script for 2D Boussinesq equations under streamfunction formulation with double-periodic and Fourier+Chebyshev grid. I am doing this in order to test the effects of rigid boundaries on the  flow of my interest.

While in both the formulations I was able to build the solver and run it successfully, the results don't seem to match. I know that the double Fourier grid results are correct but when I run same model equations with the Chebyshev formulation, although running, the model doesn't evolve in the same way.

Here is a code snippet:

1. Double Fourier - Running successful ( attached animation of how the flow evolves)

problem=de.IVP(domain,variables=['psi','theta','u','w'])

problem.substitutions['L(A)'] = 'd(A,x=2) + d(A,z=2)'  # Laplacian term
problem.substitutions['H(A)'] = 'd(A,x=4) + d(A,z=4)'  # Hyperdiffusion
problem.substitutions['J(A,B)'] = 'dx(A)*dz(B)-dx(B)*dz(A)' # Jacobian
problem.add_equation("dt(L(psi)) + dx(theta)-nu*H(psi) = - J(psi,L(psi))",condition="(nx !=0) or (nz!=0)")
problem.add_equation("dt(theta) = N2*dx(psi)-J(psi,theta)",condition="(nx !=0) or (nz!=0)")
problem.add_equation("u + dz(psi) = 0",condition="(nx !=0) or (nz!=0)")
problem.add_equation("w - dx(psi) = 0",condition="(nx !=0) or (nz!=0)")
problem.add_equation("psi=0",condition="(nx==0) and (nz==0)")
problem.add_equation("u = 0", condition="(nx==0) and (nz==0)")
problem.add_equation("w = 0", condition="(nx==0) and (nz==0)")
problem.add_equation("theta = 0", condition="(nx==0) and (nz==0)")

 2. Fourier +  Chebshev ( Running, but results not as expected)

problem=de.IVP(domain,variables=['psi','u','theta','Dtheta','zeta','Dzeta'])
problem.meta['psi','theta','zeta']['z']['dirichlet'] = True

problem.substitutions['w'] = "dx(psi)"

problem.add_equation("zeta +dx(w)-dz(u) = 0")
problem.add_equation("dt(zeta) - nu*(dx(dx(zeta)) + dz(Dzeta)) + dx(theta) = -w*Dzeta-u*dx(zeta)")
problem.add_equation("dt(theta) - (dx(dx(theta)) + dz(Dtheta)) - N2*w = -(w*Dtheta +u*dx(theta))")
problem.add_equation("u + dz(psi) = 0")
problem.add_equation("Dtheta - dz(theta) = 0")
problem.add_equation("Dzeta - dz(zeta) = 0")
problem.add_bc("left(theta) = 0")
problem.add_bc("right(theta) = 0")
problem.add_bc("left(psi) = 0")
problem.add_bc("right(psi) = 0")
problem.add_bc("left(zeta) = 0")
problem.add_bc("right(zeta) = 0")

Perhaps I am missing something very obvious here. Any suggestion is highly appreciated

[Keaton Burns]

Hi Kesava,

Here are a few things I noticed that I think all added up to the differences. First, the double-Fourier code can be optimized in a few ways:

1) Only psi and theta need to be defined as variables, and u and w can be defined using substitutions instead of additional variables and equations. This should reduce the memory footprint and linear solver time.

2) The Jacobian substitution can be replaced with an Advection substitution that references the u and w substitutions. This will reduce redundant computaitons of the derivatives of psi on the RHS and might help with performance a bit.

3) For the hyperdiffusion you probably want L(L(A)) — just doing the 4th x and z derivatives misses some cross terms.

4) Stricly speaking, a separate gauge equation only needs to be defined for psi and not theta, since psi appears in the Laplacian in the time derivative, while theta appears directly.

That leaves us with these equations:

problem = de.IVP(domain,variables=['psi','theta'])

problem.substitutions['u'] = "-dz(psi)"

problem.substitutions['w'] = "dx(psi)"

problem.substitutions['L(A)'] = 'd(A,x=2) + d(A,z=2)'  # Laplacian term

problem.substitutions['H(A)'] = 'L(L(A))'  # Hyperdiffusion term

problem.substitutions['Adv(A)'] = "u*dx(A) + w*dz(A)"

problem.add_equation("dt(L(psi)) + dx(theta) - nu*H(psi) = - Adv(L(psi))", condition="(nx!=0) or (nz!=0)")

problem.add_equation("dt(theta) = N2*w - Adv(theta)")

problem.add_equation("psi = 0", condition="(nx==0) and (nz==0)")

Now for the Fourier-Chebyshev formulation:

1) The N2 profile was different — I changed it back to match the double-Fourier version.

2) The definition of zeta here is different by a sign than in the double-Fourier code. There it is zeta = lap(psi), while here it was zeta = - lap(psi). This means the dx(theta) term has a different sign from before, so I changed the definition here to match the double-Fourier version.

3) In the double-Fourier code, theta has no diffusion, while here it has diffusion with constant 1. I think this changes the solution a lot — I instead put in nu as the diffusion constant (so Prandtl=1) and that seems to help a lot.

4) The psi initial condition was turned off here. I added that back in. It’s also important to differentiate it up to set the correct initial condition for zeta, since it appears in a time derivative in the equations.

All together the Fourier-Chebyshev equations are then:

problem=de.IVP(domain,variables=['psi','u','theta','Dtheta','zeta','Dzeta'])

problem.substitutions['w'] = "dx(psi)"

problem.substitutions['L(A,Az)'] = "dx(dx(A)) + dz(Az)"

problem.substitutions['Adv(A,Az)'] = "u*dx(A) + w*Az"

problem.add_equation("-zeta + dx(w) - dz(u) = 0")

problem.add_equation("dt(zeta) - nu*L(zeta,Dzeta) + dx(theta) = - Adv(zeta,Dzeta)")

problem.add_equation("dt(theta) - nu*L(theta,Dtheta) - N2*w = -Adv(theta,Dtheta)")

problem.add_equation("u + dz(psi) = 0")

problem.add_equation("Dtheta - dz(theta) = 0")

problem.add_equation("Dzeta - dz(zeta) = 0")

problem.add_bc("left(theta) = 0")

problem.add_bc("right(theta) = 0")

problem.add_bc("left(psi) = 0")

problem.add_bc("right(psi) = 0")

problem.add_bc("left(zeta) = 0")

problem.add_bc("right(zeta) = 0")

With all those changes, the evolution seems to be pretty similiar now.

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/751876ff-69d7-4083-bca5-c8e90f4abfecn%40googlegroups.com.

[Kesava Prasad]

Hi Keaton,

thanks for the detailed response . I really appreciate it. I apologize for commenting out certain stuff, that must have been confusing. Fascinating how even after modeling multiple problems in dedalus, you still discover something new :)

4) The psi initial condition was turned off here. I added that back in. It’s also important to differentiate it up to set the correct initial condition for zeta, since it appears in a time derivative in the equations.

I was looking through the docs here: https://dedalus-project.readthedocs.io/en/latest/notebooks/dedalus_tutorial_fields_operators.html, but couldn't find anything on why the differentiation is done in the coefficient space and not in the grid space. Is there any particular reason for it?

Thanks,

Kesava

[Daniel Lecoanet]

The derivative of a Fourier or Chebyshev series can be evaluated exactly using the Fourier/Chebyshev coefficients. If you took derivatives on the grid, you would have to use something like a finite difference approximation to the derivative, which is not exact.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/089204b3-e8f9-4bbb-b8b7-2a6b6d19753fn%40googlegroups.com.

[Keaton Burns]

To further clarify, yes the derivatives are always computed internally in coefficient space. But in the script I differentiated up the psi initial conditions like so:

u['c'] = -psi.differentiate('z')['c']

Breaking this down, "psi.differentiate('z’)” computes the derivative (spectrally), and returns a new field. I’m then copying the coefficient-space data from that field into the “u” field from the state vector, setting the initial condition for “u” in the solver.  You could also copy the data in grid space, but you might take an extra step to make sure the scales are the same between the “u” field and the new field from the derivative operator.  Copying in coefficent space just avoids having to worry about this.

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/1BEF548C-5049-4C9B-BCA5-516C8E016095%40northwestern.edu.

[Jeffrey S. Oishi]

Hi Kesava,

On Wed, Feb 2, 2022 at 10:30 AM Kesava Prasad <kesa...@gmail.com> wrote:

I was looking through the docs here: https://dedalus-project.readthedocs.io/en/latest/notebooks/dedalus_tutorial_fields_operators.html, but couldn't find anything on why the differentiation is done in the coefficient space and not in the grid space. Is there any particular reason for it?

 This is the essence of the spectral method: because we expand things on a vector space of polynomial basis functions with the property that their derivatives project (sparsely!) onto that same space, doing derivatives in coefficient space makes them essentially as accurate as possible for a given number of expansion coefficients. I suggest you have a look at section II.A.1 of the Dedalus methods paper,

https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.2.023068

And the textbook by Boyd referenced therein.

j

 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/089204b3-e8f9-4bbb-b8b7-2a6b6d19753fn%40googlegroups.com.

[Kesava Prasad]

I really appreciate the explanations. Thanks a lot.

Cheers,

Kesava