error: factor is exactly singular

[Guangyao Cui]

Hi all,

I am working on a 2D flow with external forcing in a doubly periodic domain (I attached a reference, please see equation 11-13), and I tried to add the equation in vorticity-stream function form, but it always gives an error "factor is exactly singular".

I suspect that the error comes from the definition of the equation (or condition) I defined, but I couldn't figure it out where exactly the error is coming from. I attached the script and the error message also.

Any advice is greatly appreciated. Thanks.

Guangyao

[Guangyao Cui]

I also tried to define the equation like this:

problem.add_equation("dt(om) + k**2/Re*om - 1/Re * L(om) = -J(om, psi) - Fy", condition="(nx!=0) or (ny!=0)")
problem.add_equation("L(psi) + om = 0", condition="(nx!=0) or (ny!=0)")
problem.add_equation("u - dy(psi) = 0")
problem.add_equation("v + dx(psi) = 0")
problem.add_equation("integ(integ(u,'x'), 'y')=0", condition="(nx==0) and (ny==0)")
problem.add_equation("integ(integ(v,'x'), 'y')=Uy", condition="(nx==0) and (ny==0)")

Still it doesn't work, and it shows the same error.

Guangyao

[Daniel Lecoanet]

Hi Guangyao,

You’re close! Since you’re using a Fourier series expansion in both x and y, one of the modes in the expansion varies like exp(i*0*x + i*0*y). We call that mode nx=0 and ny=0. For that mode, the x and y derivatives are both zero. In Dedalus, the left-hand side of the equation is treated implicitly, which involves a matrix inverse. Singular matrix errors occur when the matrix associated with the left-hand side terms is singular.

What are the left-hand side of the equations for the nx=ny=0 mode?

dt(om) + k^2/Re*om = 0

om = 0

u = 0

v = 0

The problem is that you have two equations for om, and no equations for psi. You’ve put conditions on some of your equations, so taking into account the conditions, we have

u = 0

v = 0

u = 0

v = 0

Now you have two equations for u and v, but no equations for om or psi.

The way to fix the issue is to put a condition on the L(psi)+om=0 equation, and replace that equation with psi=0.

Daniel

--
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/83b7bfb1-b149-4ad5-94ab-b4f08b11ebb4n%40googlegroups.com.

[Guangyao Cui]

Hi Daniel,

Thanks a lot for your response! Just one simple question:

for equations like this:

problem.add_equation("u- dy(psi) = 0")

problem.add_equation("v+ dx(psi)= 0")

Is that correct that there is no point to put a condition for u at nx=0 and ny=0, because as you explained, this equation "problem.add_equation("u- dy(psi) = 0")" will include (u=0, condition="(nx==0) and (ny==0)") by default (due to the fact that any derivative should be zero at nx=0, ny=0)? Therefore I just need to put conditions for om, and psi. Is that correct? 

Thanks,

Guangyao

You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/94mqhXOXyOs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/D34A2AFE-E631-4B52-B613-8583168D7E2C%40northwestern.edu.

[Daniel Lecoanet]

Yes, that’s right. And you don’t need a condition on the dt(om) + … equation because that one still makes sense when nx=ny=0.

Daniel

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

[Guangyao Cui]

I see. Thanks a lot.

Guangyao

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/25CC69D1-54E3-4938-96C6-C8142C93C853%40northwestern.edu.

[Eric Tseng]

Hi all, 

I was working on a similar problem but on a spherical coordinate. 

Part of my code looks like this:

# Parameters

Nphi = 256

Ntheta = 128

dealias = 3/2

R   = 6.37122e6

coords = de.S2Coordinates('phi', 'theta')

dist      = de.Distributor(coords, dtype=dtype)

basis   = de.SphereBasis(coords, (Nphi, Ntheta), radius=R, dealias=dealias, dtype=dtype)

...

problem.add_equation("lap(psi_b) = vor_b" , condition='(Ntheta!=0)')

problem.add_equation("psi_b = 0" , condition='(Ntheta==0)')

However, I got the following error message when I tried to solve the inverse problem 

NameError: name 'Ntheta' is not defined

I thought I have defined Ntheta at the beginning?

Any suggestion is appreciated. 

Thx!

Eric  

[Keaton Burns]

Hi Eric,

Here “Ntheta” is the resolution of the basis. Internally “n{coord_name}” is used for equation condition selections, so you should use “ntheta”. However, we recommend using scalar variables to help enforce gauge conditions, rather than equation conditions, which you can read more about here: https://dedalus-project.readthedocs.io/en/latest/pages/gauge_conditions.html 

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/29f8a520-6f8e-4fa9-8198-f351d80ddaa2n%40googlegroups.com.