Singular matrix error after two iterations and overflow error

[jc2...@cornell.edu]

Hi,

I am working on a 2D NS equation in streamfunction-vorticity formulation with topography in a tunnel (Fourier-Chebyshev base). I have built a similar case in double periodic and it seem to work. However, when migrating the case to the new mixed bases with some revised boundary conditions, the main loop does 2 iterations, then report "overflow encountered in multiply" warning and then "factor is exactly singular" error in the next iteration. I have encountered the singular matrix error previously, and it is usually related to improper boundary condition. However, I am not sure what is the probable cause for the singular matrix error after several iteration. Is it still the B.C. problem, or is it because of the overflow and thus make the matrix singular?

My governing equation and boundary conditions are as follow:

problem = de.IVP(domain, variables=['Psi', 'zeta','Psix','Psiy'])

problem.parameters['nu'] = nu 

problem.parameters['H'] = H

problem.parameters['Hinv'] = Hinv

problem.parameters['r'] = r

problem.substitutions['v'] = "Hinv*dx(Psi)"

problem.substitutions['u'] = "-Hinv*dy(Psi)"

problem.add_equation("Psix-dx(Psi)=0")

problem.add_equation("Psiy-dy(Psi)=0")

problem.add_equation("dt(zeta) + r*(Hinv**2)*(dx(Psix)+dy(Psiy))= Psiy*dx(zeta*Hinv) - Psix*dy(zeta*Hinv) + nu*(dx(dx(zeta)) + dy(dy(zeta)))+r*(Hinv**3)*(Psix*dx(H))+F()")

problem.add_equation("zeta  = dx(v) - dy(u)")

problem.add_bc("left(Psi)=12")

problem.add_bc("right(Psi)=12")

Any suggestion or advice on how to fix this is highly appreciated!

Best,

Jinshi

[Evan Henry Anders]

Hi Jinshi,

I can't say for sure what the problem is, but it sounds like something about the equations as written are numerically unstable, and that instability is showing up over the course of a few timesteps. When it shows up, a multiply operation is overflowing, making some data NaN or inf which then leads to a singular matrix on the next timestep I think. This behavior makes me think it's something going wrong with the right-hand side of the equations (the explicit terms), which are more sensitive to timestep choices, etc.

One quick thing I'm seeing that's almost definitely going wrong: diffusion terms should be timestepped implicitly, not explicitly (diffusion terms are unstable for explicit timestepping methods). This means that you want to move  the "nu*(dx(dx(zeta)) + dy(dy(zeta)))" term to the left-hand side of the equations so that it timesteps implicitly. Since you're now trying to make this work in a Fourier-Chebyshev domain, this probably requires an equation that defines "zeta_y - dy(zeta) = 0," and this inclusion on the LHS likely also requires two additional boundary conditions.

Hopefully this is helpful as a starting point. Best,

Evan

--
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/81433799-4b2a-4263-962e-675c5726a9e1n%40googlegroups.com.

[Jeffrey S. Oishi]

Just a quick note:

On Tue, Feb 9, 2021 at 7:28 PM Evan Henry Anders <Evan....@colorado.edu> wrote:

One quick thing I'm seeing that's almost definitely going wrong: diffusion terms should be timestepped implicitly, not explicitly (diffusion terms are unstable for explicit timestepping methods). 

It is certainly true that you probably want to have your diffusion terms on the left, not right side as Evan says. However, it is not true that diffusion terms are unstable for explicit timestepping methods. They are just much stiffer, requiring much smaller timesteps.

 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CALW8q5OCXMAFq0uON2egyd9jvaddzJC4_wv5cD%2BPSdEirFJizw%40mail.gmail.com.

[jc2...@cornell.edu]

Hi Evan and Jeff,

Thank you so much for your suggestions! I have managed to get it work after moving the diffusion term to the right with proper extra boundary conditions. 

Best,

Jinshi