Re: matrix is singular for 1D compressible NS

[Daniel Lecoanet]

Hi,

I would suggest that you look at an example of how we've implemented compressible hydrodynamics in Dedalus in the past, e.g., http://adsabs.harvard.edu/abs/2014ApJ...797...94L.  I don't think your formulation of the equations will work in Dedalus because of the way we do implicit-explicit timestepping and apply boundary conditions.

Daniel

On Tue, Nov 28, 2017 at 1:21 PM, Xin Bian <bianx...@gmail.com> wrote:

Hi Dedalus Group,

I just start using Dedalus.

I started with 2D compressible NS, but it says matrix is singular. Then I tried to do it in 1D but got the same problem.

Another thing is, I intended to implement dz(rho) = 0 on both bottom and top, w=0 on both bottom and top. It seems that I should specify 3 conditions for three pdes.

The form of compressible NS and the code are attached. 

In the code, heat flux in energy equation is written in this form:

Thanks.

Xin

--
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-users+unsubscribe@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/548266d5-e0d7-4925-9d14-31a3d440d05d%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Xin Bian]

Thanks Daniel. I'll look at that

On Tuesday, November 28, 2017 at 1:49:37 PM UTC-5, Daniel Lecoanet wrote:

Hi,

I would suggest that you look at an example of how we've implemented compressible hydrodynamics in Dedalus in the past, e.g., http://adsabs.harvard.edu/abs/2014ApJ...797...94L.  I don't think your formulation of the equations will work in Dedalus because of the way we do implicit-explicit timestepping and apply boundary conditions.

Daniel

On Tue, Nov 28, 2017 at 1:21 PM, Xin Bian <bianx...@gmail.com> wrote:

Hi Dedalus Group,

I just start using Dedalus.

I started with 2D compressible NS, but it says matrix is singular. Then I tried to do it in 1D but got the same problem.

Another thing is, I intended to implement dz(rho) = 0 on both bottom and top, w=0 on both bottom and top. It seems that I should specify 3 conditions for three pdes.

The form of compressible NS and the code are attached. 

In the code, heat flux in energy equation is written in this form:

Thanks.

Xin

--
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.

[Xin Bian]

Hi Daniel,

I looked at that your paper. You decomposed the field into background and fluctuation. I am not sure if I understand the paper correctly. Is that because of the heat conduction acting different on perturbation fields and background field or the requirement of the code?

If it is because of the physics in heat conduction, I don't care about that part for now. I put 2D version of the following equations in Dedalus and assume constant diffusivity:

 

from  Anders, Evan H., and Benjamin P. Brown. "Convective heat transport in stratified atmospheres at low and high Mach number." Physical Review Fluids 2.8 (2017): 083501.

Then I specified w = dz(T) = u = 0 at top and bottom as BCs. I also want to specify dz(rho) = 0, but not sure how to do that. For not, I just make rho = 5 at top and rho = 2 at bottom (the same as ICs).

The matrix can be generated, but the code blows up immediately. This is the IC at t = 0.

 

The example is attached. Did I do something wrong here?

Thanks.

Xin

On Tuesday, November 28, 2017 at 1:49:37 PM UTC-5, Daniel Lecoanet wrote:

Hi,

I would suggest that you look at an example of how we've implemented compressible hydrodynamics in Dedalus in the past, e.g., http://adsabs.harvard.edu/abs/2014ApJ...797...94L.  I don't think your formulation of the equations will work in Dedalus because of the way we do implicit-explicit timestepping and apply boundary conditions.

Daniel

On Tue, Nov 28, 2017 at 1:21 PM, Xin Bian <bianx...@gmail.com> wrote:

Hi Dedalus Group,

I just start using Dedalus.

I started with 2D compressible NS, but it says matrix is singular. Then I tried to do it in 1D but got the same problem.

Another thing is, I intended to implement dz(rho) = 0 on both bottom and top, w=0 on both bottom and top. It seems that I should specify 3 conditions for three pdes.

The form of compressible NS and the code are attached. 

In the code, heat flux in energy equation is written in this form:

Thanks.

Xin

--
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.

[Evan H. Anders]

Hey Xin,

Usually we solve on these 7 variables in the 2D FC equations:

['u','u_z','w','w_z','T1', 'T1_z', 'ln_rho1']

Then using substitutions we specify the different components of the stress tensor, e.g.,

problem.substitutions["txx"] = "(2*dx(u) - 2/3*Div_u)"

problem.substitutions["tyy"] = "(2*dy(v) - 2/3*Div_u)"

problem.substitutions["tzz"] = "(2*w_z   - 2/3*Div_u)"

which lets us write the momentum equation in a similar form to how you wrote it.  This system has 7 equations:

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

problem.add_equation("dz(w) - w_z = 0")

problem.add_equation("dz(T1) - T1_z = 0")

+ momentum, x

+ momentum, z

+ energy

+ continuity

and then our boundary conditions are: 2x thermal, 4x velocity (impenetrable, and then no-slip or stress-free).

The system you're solving right now has 8 equations, but only 7 unique z-derivatives in it.  In general, a rule of thumb you have to follow in Dedalus is "max number of boundary conditions I can use = number of unique chebyshev (z) derivatives I have in my system."  So here in your system, the 8th boundary condition over-constrains it and you get a singular matrix.  It spits this out as a warning when you I it, but it doesn't outright fail, and that's why it explodes immediately.

As far as setting the density constant at the boundaries, I think you're going to have to set your velocity boundary conditions so that d ln rho / dt = 0  by definition at the boundaries (so plug your velocity BCs into the continuity equation and see if that enforces your constraint of not letting rho change at the boundaries).

Hope this helps!

Evan

To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-users+unsubscribe@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/079d62ad-6b07-4b76-ba0a-03301458feb4%40googlegroups.com.

[Xin Bian]

Thanks Evan.

I should use substitution rather than equation for the stress tensor.

Now I have these 8 variables, I added extra one rho_1 as density for the use in BCs.

problem = de.IVP(domain, variables=['w', 'w_z','u', 'u_z', 'log_rho','T', 'T_z', 'rho_1'])

problem.substitutions["t_xx"] = ("dim_cof*(2*dx(u) - w_z)")

problem.substitutions["t_xz"] =("dx(w) + u_z")

problem.substitutions["t_zz"] =("dim_cof*(2*w_z - dx(u))")

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

problem.add_equation("dz(w) - w_z = 0")

problem.add_equation("dz(T) - T_z = 0")

problem.add_equation("rho_1 = exp(log_rho)")

problem.add_equation("dt(log_rho) + w_z + dx(u) = -w*dz(log_rho) - u*dx(log_rho)")

problem.add_equation("dt(w) + T_z   - nu*dz(t_zz) - nu*dx(t_xz) = - w*w_z - u*dx(w) - T*dz(log_rho) + g + nu*t_zz*dz(log_rho) + nu*t_xz*dx(log_rho)")

problem.add_equation("dt(u) + dx(T) - nu*dx(t_xx) - nu*dz(t_xz) = - w*dz(u) - u*dx(u) - T*dx(log_rho)   + nu*t_xx*dx(log_rho) + nu*t_xz*dz(log_rho)")

problem.add_equation("dt(T) - 1/Cv*(kappa*(dz(T_z) + dx(dx(T)))) = - w*T_z - u*dx(T) - (gamm-1)*T*(w_z + dx(u)) + 1/Cv*(kappa*(T_z*dz(log_rho) + dx(T)*dx(log_rho)) + nu*t_zz*w_z + nu*t_xx*dx(u) + nu*t_xz*u_z + nu*t_xz*dx(w))")

And BCs:

problem.add_bc("right(dz(rho_1)) = 0")

problem.add_bc("left(dz(rho_1)) = 0")

problem.add_bc("left(dz(u)) = 0")

problem.add_bc("right(dz(u)) = 0")

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

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

I tried a simple hydrostatic IC with constant density=4 and it looks like:

I still got NAN after 1st or several steps. Then I tried several combination of BCs, like u=w=T_z=0, u=w=dz(rho_1)=0. 

The only combination that doesn't blow up is g = 0 and u=w=T_z=0, where everything is just constant.

To locate the problem, I also tried Euler equations by setting kappa = nu = 0, still got the same thing.

Are there other things I should consider? If not, maybe I have some bugs in the code.

Thanks.

[Daniel Lecoanet]

Hi Xin,

It's probably a good idea to start off with a set of equations that work, such as the equations in Evan's paper or my paper.  Can you get those equations to work?

Daniel

On Wed, Nov 29, 2017 at 5:21 PM, Xin Bian <bianx...@gmail.com> wrote:

Thanks Evan.

I should use substitution rather than equation for the stress tensor.

Now I have these 8 variables, I added extra one rho_1 as density for the use in BCs.

problem = de.IVP(domain, variables=['w', 'w_z','u', 'u_z', 'log_rho','T', 'T_z', 'rho_1'])

problem.substitutions["t_xx"] = ("dim_cof*(2*dx(u) - w_z)")

problem.substitutions["t_xz"] =("dx(w) + u_z")

problem.substitutions["t_zz"] =("dim_cof*(2*w_z - dx(u))")

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

problem.add_equation("dz(w) - w_z = 0")

problem.add_equation("dz(T) - T_z = 0")

problem.add_equation("rho_1 = exp(log_rho)")

problem.add_equation("dt(log_rho) + w_z + dx(u) = -w*dz(log_rho) - u*dx(log_rho)")

problem.add_equation("dt(w) + T_z   - nu*dz(t_zz) - nu*dx(t_xz) = - w*w_z - u*dx(w) - T*dz(log_rho) + g + nu*t_zz*dz(log_rho) + nu*t_xz*dx(log_rho)")

problem.add_equation("dt(u) + dx(T) - nu*dx(t_xx) - nu*dz(t_xz) = - w*dz(u) - u*dx(u) - T*dx(log_rho)   + nu*t_xx*dx(log_rho) + nu*t_xz*dz(log_rho)")

problem.add_equation("dt(T) - 1/Cv*(kappa*(dz(T_z) + dx(dx(T)))) = - w*T_z - u*dx(T) - (gamm-1)*T*(w_z + dx(u)) + 1/Cv*(kappa*(T_z*dz(log_rho) + dx(T)*dx(log_rho)) + nu*t_zz*w_z + nu*t_xx*dx(u) + nu*t_xz*u_z + nu*t_xz*dx(w))")

And BCs:

problem.add_bc("right(dz(rho_1)) = 0")

problem.add_bc("left(dz(rho_1)) = 0")

problem.add_bc("left(dz(u)) = 0")

problem.add_bc("right(dz(u)) = 0")

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

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

I tried a simple hydrostatic IC with constant density=4 and it looks like:

I still got NAN after 1st or several steps. Then I tried several combination of BCs, like u=w=T_z=0, u=w=dz(rho_1)=0. 

The only combination that doesn't blow up is g = 0 and u=w=T_z=0, where everything is just constant.

To locate the problem, I also tried Euler equations by setting kappa = nu = 0, still got the same thing.

Are there other things I should consider? If not, maybe I have some bugs in the code.

Thanks.

On Wednesday, November 29, 2017 at 1:53:44 PM UTC-5, Evan H. Anders wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/079d62ad-6b07-4b76-ba0a-03301458feb4%40googlegroups.com.

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

--
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-users+unsubscribe@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/cc4a68a9-3900-41d3-be8f-ad844e5adbde%40googlegroups.com.

[Xin Bian]

Hi Daniel,

I'll try that. Thanks.