Matrix is exactly singular

[Jared Whitehead]

I'm running into this error, and I'm sure it is because I am setting things up incorrectly.  I'm trying to simulate 2D Boussinesq (NS + stably stratifying active scalar) in a doubly periodic domain (I'm attaching the modified script I generated from the 2d RB example).

My first guess is that I'm missing something for the doubly periodic boundary conditions...is it enough to just specify the domain is Fourier?

[Daniel Lecoanet]

problem.add_equation("thetaz - dx(theta) = 0")

Should be 

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

Otherwise looks like it should work.

On Fri, May 29, 2015 at 6:18 PM, Jared Whitehead <whit...@mathematics.byu.edu> wrote:

I'm running into this error, and I'm sure it is because I am setting things up incorrectly.  I'm trying to simulate 2D Boussinesq (NS + stably stratifying active scalar) in a doubly periodic domain (I'm attaching the modified script I generated from the 2d RB example).

My first guess is that I'm missing something for the doubly periodic boundary conditions...is it enough to just specify the domain is Fourier?

--
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 post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/b7f55e2a-7221-4430-8bb4-0e2a7dafcbff%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jared Whitehead]

Daniel-

Thanks for catching that line.  Unfortunately I still get the same error even with the correction.  Any ideas?

Jared Whitehead

Assistant Professor, Mathematics Department
Brigham Young University

350 TMCB, Provo, UT, USA

http://math.byu.edu/~whitehead

http://www.mormon.org/me/15TK

--
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/Syupo6xWvsU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.

To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAJoYf%3DhfhKJhwg-TnRkMJnHEYOdAkyY3%3Dtr7MrTd6Cz7Jz2U1g%40mail.gmail.com.

[Ben Brown]

Jared,

      First, you need to remove all of the "psix - dx(psi) = 0" eqations: the first order reduction is done only in the Chebyshev direction.  That's done in the attached 2d_bouss_singular.py

This still fails with a

ValueError: Pencil [0, 0] has 6 equations for 7 variables.

      This second problem is, I think, evidence that you're running into a classic problem with constrained equations (e.g., equation sets with div.u = 0).  Even though this is doubly periodic and there are no boundary conditions here, one of the equations becomes redundant for the nx=0 mode, where the divergence constraint leads to dz(psi) = 0 etc.

I did a version of these 2-D periodic equations for boussinesq KH here:

      https://bitbucket.org/bpbrown/kelvin_helmholtz/src/tip?at=default

This was with the old parser, and there you can see that, at line 43 of equations.py, I have

self.problem.add_int_bc("P = 0", condition="dx == 0")

Here we set the gauge of the lagrange multiplier (P) to be zero as our non-redundant equation.

We need the equivalent line for your problem, for the new parser, where now we don't have add_int_bc, and dx -> nx.

I'm unclear right now what the equivalent lagrange multiplier is in this streamfunction system.  I tried adding the following:

problem.add_equation("omega = 0", condition="nx==0")

and some variants, but am still getting a mismatch between N_equations and N_variables for some of the low pencils ([0,0], [1,0], [0,1] or [1,1] variously).

Anyone have further input on this?  I'll keep playing with it some later this morning.

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAFHSqE%3DnrGhWv%2BO1BHX%2BNtLfF4KZ%2B7w3PgqUk2obx6VTjXjGsw%40mail.gmail.com.

[Keaton Burns]

Hi Jared,

For the nx=nz=0 pencil, dx=dz=0 and so the first and second equations end up being redundant constraints on omega.  To fix this, you can replace the first equation with the following two conditional equations, which instead fixes the gauge of the stream function for the constant pencil:

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

problem.add_equation("omega - dx(dx(psi)) - dz(dz(psi)) = 0", condition="(nx!=0) or (nz!=0)")

As Ben noted, first-order reductions aren’t necessary for Fourier derivatives, but they also shouldn’t prevent the code from running.  Since you’re doing doubly periodic, I’ve removed both sets of derivative reductions in the version below.  I’ve also thrown in some parsing substitutions defining “u” and “w” in terms of the stream function, which are completely optional but can be used in the equations, analysis tasks, etc.

Finally, there was also a slight bug in the Fourier matrix constructor which I recently introduced, but I’ve fixed it, so the attached version should work with an updated copy of Dedalus.  

-Keaton

[Jared Whitehead]

Thanks for the help guys, and sorry it is taking me so long to try and implement these changes.

Keaton, this is what happens when I use your script.  I tried the RK222 solver as well and got the same error...

Is there something missing in the way I set up the domain?

(dedalus)new-host-3:Downloads knrider$ python3 2d_bouss.py 

2015-06-02 21:53:51,905 pencil 0/1 INFO :: Building pencil matrix 1/32 (~3%) Elapsed: 0s, Remaining: 1s, Rate: 2.1e+01/s

2015-06-02 21:53:52,051 pencil 0/1 INFO :: Building pencil matrix 4/32 (~12%) Elapsed: 0s, Remaining: 1s, Rate: 2.1e+01/s

2015-06-02 21:53:52,248 pencil 0/1 INFO :: Building pencil matrix 8/32 (~25%) Elapsed: 0s, Remaining: 1s, Rate: 2.1e+01/s

2015-06-02 21:53:52,444 pencil 0/1 INFO :: Building pencil matrix 12/32 (~38%) Elapsed: 1s, Remaining: 1s, Rate: 2.0e+01/s

2015-06-02 21:53:52,638 pencil 0/1 INFO :: Building pencil matrix 16/32 (~50%) Elapsed: 1s, Remaining: 1s, Rate: 2.1e+01/s

2015-06-02 21:53:52,833 pencil 0/1 INFO :: Building pencil matrix 20/32 (~62%) Elapsed: 1s, Remaining: 1s, Rate: 2.1e+01/s

2015-06-02 21:53:53,027 pencil 0/1 INFO :: Building pencil matrix 24/32 (~75%) Elapsed: 1s, Remaining: 0s, Rate: 2.1e+01/s

2015-06-02 21:53:53,223 pencil 0/1 INFO :: Building pencil matrix 28/32 (~88%) Elapsed: 1s, Remaining: 0s, Rate: 2.1e+01/s

2015-06-02 21:53:53,419 pencil 0/1 INFO :: Building pencil matrix 32/32 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 2.0e+01/s

2015-06-02 21:53:53,421 __main__ 0/1 INFO :: Solver built

2015-06-02 21:53:53,422 __main__ 0/1 INFO :: Starting loop

2015-06-02 21:53:53,505 __main__ 0/1 ERROR :: Exception raised, triggering end of main loop.

2015-06-02 21:53:53,505 __main__ 0/1 INFO :: Iterations: 0

2015-06-02 21:53:53,505 __main__ 0/1 INFO :: Sim end time: 0.000000

2015-06-02 21:53:53,505 __main__ 0/1 INFO :: Run time: 0.08 sec

2015-06-02 21:53:53,505 __main__ 0/1 INFO :: Run time: 0.000023 cpu-hr

Traceback (most recent call last):

  File "2d_bouss.py", line 96, in <module>

    solver.step(dt)

  File "/Users/knrider/Desktop/dedalus/src/dedalus/dedalus/core/solvers.py", line 339, in step

    self.timestepper.step(self, dt, wall_time)

  File "/Users/knrider/Desktop/dedalus/src/dedalus/dedalus/core/timesteppers.py", line 128, in step

    MX0.set_pencil(p, p.M*x)

AttributeError: 'Pencil' object has no attribute 'M'

Jared Whitehead

Assistant Professor, Mathematics Department
Brigham Young University

350 TMCB, Provo, UT, USA

http://math.byu.edu/~whitehead

http://www.mormon.org/me/15TK

--
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/Syupo6xWvsU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/C6A53310-6E0A-42C9-BB0B-726A98AA0010%40gmail.com.

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

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

<2d_bouss_singular.py>

--
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/Syupo6xWvsU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/C6A53310-6E0A-42C9-BB0B-726A98AA0010%40gmail.com.

[Keaton Burns]

Hi Jared,

Nope, this is my mistake -- another Fourier-only matrix construction bug slipped by me. It should work now.

Sorry!
-Keaton

> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAFHSqEmvFnfwJozZK2b%2BsAow7WC_7B1oHn1U%3DzEw2R1wVPdNyg%40mail.gmail.com.

[Jared Whitehead]

Alright, looks like it ran, and the output looks like I would expect.  Now I'll see if I can mess things up again :)

Thanks!

Jared Whitehead

Assistant Professor, Mathematics Department
Brigham Young University

350 TMCB, Provo, UT, USA

http://math.byu.edu/~whitehead

http://www.mormon.org/me/15TK

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/689D08C7-3009-4AF2-82B5-A257B5D5B8CE%40gmail.com.