Singular matrix misery

[Daniel Lecoanet]

Hi all,

I'm getting singular matrices and don't know why.

I've attached a script where I have the following problem:

linear_wave_problem.add_equation("dt(u) + dx(pomega) - nu*( dx(dx(u)) + dy(dy(u)) + dz(u_z) ) = - u*dx(u) - v*dy(u)")

linear_wave_problem.add_equation("dt(v) + dy(pomega) - nu*( dx(dx(v)) + dy(dy(v)) + dz(v_z) ) = - u*dx(v) - v*dy(v)")

linear_wave_problem.add_equation("dx(u) + dy(v) = 0", condition = "dx != 0 or dy != 0")

linear_wave_problem.add_equation("pomega = 0", condition = "dx == 0 and dy == 0")

linear_wave_problem.add_equation("dz(u) - u_z = 0")

linear_wave_problem.add_equation("dz(v) - v_z = 0")

linear_wave_problem.add_left_bc( "u = 0")

linear_wave_problem.add_right_bc("u = 0")

linear_wave_problem.add_left_bc( "v = 0")

linear_wave_problem.add_right_bc("v = 0")

All the LHS matrices are singular.  If you look at the LHS matrices (solver.pencils[i].LHS.todense()), they are not obviously singular, i.e., there's an entry in every column and row.

Any thoughts?

Daniel

[Ben Brown]

Daniel,

    Why does your div(u) equation have a condition on it?

--Ben

--
You received this message because you are subscribed to the Google Groups "Dedalus Development" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-dev...@googlegroups.com.
To post to this group, send email to dedal...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgWBeB_fL1KRYFi3zzwz4yh%3DepYDwyhZt%2BC6x5qLT_%2BiA%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

[Keaton Burns]

Hi Daniel,

The LHS matrices are linear combinations of L and M based on the variable timestep, and so they’re just empty before you start time stepping.  Instead you should probably look at the L matrix — and by examining the condition numbers it looks like the dx=dy=0 matrix is not singular, but the other ones are.

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Development" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-dev...@googlegroups.com.
To post to this group, send email to dedal...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgWBeB_fL1KRYFi3zzwz4yh%3DepYDwyhZt%2BC6x5qLT_%2BiA%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

<test_problem.py>

[Daniel Lecoanet]

The reason why I looked at LHS in the first place is that timestepping is giving me singular matrices.  I've attached a new version of the script.  Now I take a timestep.  Some of the pencils have singular matrices, and others do not.  Interestingly, I found a pair of pencils that have kx and ky interchanged -- one has a singular LHS matrix, and the other does not.  I think the problem is symmetric to kx <-> ky, so I'm very confused.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/C4C8B013-F420-4DAB-8D8F-97896A600DA1%40gmail.com.

[Daniel Lecoanet]

When kx=ky=0, the div(u) equation is just 0=0 which is singular.  Since there's no dz(pomega) in this problem, you need to set a pomega gauge at every z level, i.e., you need a special equation for pomega, rather than a special boundary condition.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAHqBLzy-YKRMYqxBsouJ73vm7Y3nehsmGx%3DhbiYTvTB00a1nsw%40mail.gmail.com.

[Keaton Burns]

I think you’re just seeing an artifact of the numerical inversion.  Printing the LHS condition number for those pencils (np.linalg.cond) shows that they’re both equal and extremely large/singular (6.022e16).

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgZkQuCT7Rg%3DyCUXBrMMoz8d8_fWJvtDoAmY77wvDqtyQ%40mail.gmail.com.

[Keaton Burns]

Oh wait I had a typo — they’re not precisely equal but still extremely large (6.022e16 vs 1.14e17), so I still think they’re both singular and the difference is possibly numerical.

[Daniel Lecoanet]

Ok, I could believe that.  But I don't expect these matrices to be singular in the first place, and yet a bunch of them are!

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/499767F5-824C-4F31-B928-67C40DC3CE2A%40gmail.com.

[Daniel Lecoanet]

Oh, I should also say that if you remove the div(u) equation and remove the grad(p) from the momentum equations, then everything works properly.  So the singularness of the equations presumably is related to the div(u) equation.

Daniel

[Daniel Lecoanet]

I've altered the script to do timestepping.  When I run it on my laptop I get

0.998139819859                                                                                  

7.28447832692

97.0246069245

1261.7669741

15918.3780139

196078.2414

2369249.47983

28178815.2465

330723896.243

3837721944.12

44105960449.3

This is the amplitude of a field that's supposed to be viscously decaying.  Interestingly, if I increase the timestep, things stabilize, and the amplitude decays with time.  On the Discover computer, I get nan after the first timestep.  In my experience, different computers have different ways of dealing with singular matrices, so it's not a surprise that I get a different answer on my laptop than on Discover.

Daniel

[Ben Brown]

Daniel,

    Quick, and probably naive question:

Why does this system know anything about the z-direction? (why isn't it strictly 2-D?).

I'm not seeing anything here that will feed dynamics in the z-direction, except your initial condition, and that just couples in viscously... right?

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgW_qiujBj4pG0xCgCzF3RQ-h_%3DYzMFZPfOSwoVTSX4bQ%40mail.gmail.com.

[Daniel Lecoanet]

Viscosity acts in the z direction.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAHqBLzyj0YrZidqrP-OJAwqEdWRF2PLyv%3DgbL16HozrWHWU-rw%40mail.gmail.com.

[Daniel Lecoanet]

I wanted to make sure this problem was well-posed...  You can rewrite u and v in terms of a streamfunction, i.e.,

u = dy(psi)

v = -dx(psi)

where psi=psi(x,y,z).  Then taking the curl of the momentum equation, you get

dt (dx^2 + dy^2) psi - nu dz^2 (dx^2 + dy^2) psi = 0,

where I've dropped the x and y parts of the viscosity, and the nonlinear terms, for simplicity.  I coded up this equation and get the right answer -- viscous decay.

So I'm pretty sure the problem is theoretically well-posed, although it's certainly possible that my implementation in primitive variables is incorrect.

Daniel

[Ben Brown]

Daniel,

    What happens if you drop your pomega gauge equation from the primitive set?

--B

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DicqXPDuvSwNAGSn%3DDVbqsf8oRCEuGQy7saAeVO%2Be6%3Dcw%40mail.gmail.com.

[Daniel Lecoanet]

Same as before.  

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAHqBLzzdCL8cG4h41pQy3%2BO0ogLFg1zbRNckH%3Ds%2BmYGsh9PGmw%40mail.gmail.com.

[geoff]

I agree that your primitive system seems perfectly fine in general. I don’t think your stream-function system follows correctly from the primitive one. But even if you fix it, I think the problem should persist.  

Parsing out the primitive system, I think this is what happens:

M = 

[1 , 0 , 0 , 0 , 0]

[0 , 1 , 0 , 0 , 0]

[0 , 0 , 0 , 0 , 0]

[0 , 0 , 0 , 0 , 0]

[0 , 0 , 0 , 0 , 0]

L =

[nu*(k^2+l^2)    , 0            , i*k , -nu*dz , 0    ]

[0               , nu*(k^2+l^2) , i*l , 0      , -nu*dz]

[i*k             , i*l          , 0   , 0      , 0    ]

[ dz             , 0            , 0   , -1     , 0    ]

[ 0              , dz           , 0   , 0      , -1   ]

 

This should form the collapsed system:

det( M*dt + L ) = 

 (k^2+l^2)*{ dt + nu*(k^2+l^2 - dz^2) }

This should only complain at k=l=0. 

So, that’s all to say I don’t know what’s going on, but I checked that your not missing something simple. Unless I’m missing something simple. 

-G   

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DicqXPDuvSwNAGSn%3DDVbqsf8oRCEuGQy7saAeVO%2Be6%3Dcw%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

<test_problem_psi.py>

[Daniel Lecoanet]

On Thu, Nov 6, 2014 at 3:49 PM, geoff <geoffrey...@gmail.com> wrote:

I agree that your primitive system seems perfectly fine in general. I don’t think your stream-function system follows correctly from the primitive one. But even if you fix it, I think the problem should persist.  

Why do you think that the stream-function systems follows from the primitive one?  I neglected the diffusion in the x and y direction for simplicity because I tested that this did not solve any problems for the primitive equations.  If you do this for the primitive system, then the determinant matches the equation for psi, I think.

 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/0BF19006-6429-4417-9A77-4031A787D7C0%40gmail.com.

[geoff]

I figured that you "neglected the diffusion in the x and y direction for simplicity”, but I was commenting just in case.   The real point is that totally correct, or with a little fudge, both scaler-variable approaches should give sensible answers.  

I think I might have had a similar issue a little while ago with different MHD.

It looks like you are having problems with the algebraic constraints.  Looking at the eigenvalues, you might find part of the spectrum that is supposed to be infinity, becoming finite and large, but with the sign that would blow up the time stepping.  It all relates to zeros not actually being zeros in the generalised part of a generalised eigenvalue problem.  

I don’t know what to expect, but you might try to put time derivatives around the divergence condition. That just moves it from the L —> M operators.   Might do something? 

-G

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgMfdVcaeFnVbaVqDX7owdv-fadGoXPA2OCsQOyKaHxKw%40mail.gmail.com.

[geoff]

get on Skype if you are around. 

On 7 Nov 2014, at 1:02 pm, Daniel Lecoanet <dlec...@berkeley.edu> wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgMfdVcaeFnVbaVqDX7owdv-fadGoXPA2OCsQOyKaHxKw%40mail.gmail.com.

[Keaton Burns]

I think there’s an issue with the primitive formulation.  At least consider the case where dx = 0, dy != 0.

Then the divergence constraint becomes just “dy*v = 0”.  That’s an algebraic problem, with N equations that just depend on the N coefficients of v.  But then you also have the two boundary conditions “left(v) = 0” and “right(v) = 0”, which are two more equations that only depend on the N coefficients of v.  So the LHS matrix is singular from that: there are N+2 rows (equations) that only depend on v (N columns).

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CE834B69-D4C3-4835-A773-C31E640A4499%40gmail.com.

[Keaton Burns]

[geoff]

For posterity on the dev thread. Daniel, Keaton, and I just figured out that the system is only 2nd order in dz, but he’s got 4 boundary conditions.  

Without using a stream function (which would also work well), one could use:

dt(u) + dx(p) - nu*( dx(dx(u) + dy(dy(u)) + dz(uz) ) = F

dz(u) - uz  = 0

dx(u) + dy(v) = 0

dx(dx(p)) + dy(dy(p)) = dx(F) + dy(G) 

with boundary condition 

dx(v) - dy(u) = whatever, at z=0,1

The time-dependent v equation is normally:

dt(v) + dy(p) - nu*(dx(dx(u) + dy(dy(u)) + dz(dz(v)) = G

This gets dropped in favour of the pressure laplacian equation, which takes care of the divergence conditions without add more z derivatives. 

All of this requires dealing with dx,dy=0, cases by hand. This is the usual headache, but it’s doable.

Now we have a case of non-singular matrix misery.  

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/54100755-48A9-4C0F-8A60-C4F371AAC3AB%40gmail.com.

[Daniel Lecoanet]

Hi all,

I've attached a problem file showing an implementation of this that actually works.  The boundary conditions are a bit of a pain in the ass.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/F42EA902-D797-4758-BA2D-8E59B86D84C7%40gmail.com.