Inconsistency in the divergence

[Miguel Beneitez]

Hi,

I am working with 2D simulations of incompressible channel flow and I am trying to validate and test and test the performance of other codes against Dedalus.

To test a solution obtained with a different code, I interpolate it into the Dedalus gridspace, which causes small errors in the divergence of the field. I have computed the divergence during runtime in two different ways:

(1)  -> 

f_div_op = de.operators.differentiate(fchannel.u,'x')+de.operators.differentiate(fchannel.v,'y')

(2) ->

flow = flow_tools.GlobalFlowProperty(fchannel.solver)

flow.add_property('dx(u)+vy', name='divu_eq')

(Note that this is equivalent to the continuity equation being solved)

However both methods give different values for the divergence. It would be great if you could point out to where is the difference in the values. Thanks a lot!

Best,

Miguel

[Daniel Lecoanet]

Hi,

I assume one of your equations is something like “vy - dy(v) = 0”. Based off that, it would be very reasonable for you to assume that vy is equal to the y derivative to v. However, that is not quite the case, because you cannot simultaneously satisfy the equation and the boundary conditions. Thus vy will deviate from the y derivative of v in order to satisfy the boundary conditions. For typical simulations, the size of the discrepancy is related to the truncation error of the simulation, so it should be small. For more details, see, e.g., the section "Solving differential equations with spectral methods” on page three of the methods paper: https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.2.023068

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/6acd6b9f-397c-4075-90e0-aa40fbebf3e4n%40googlegroups.com.

[Miguel Beneitez]

Hi Daniel,

Thank you for the answer. I understand now the reason for the underlying difference!

Best,

Miguel

[Pierre-Yves Goffin]

Hi,

I had the same question as Miguel and it's also clearer for me now. But now I'm wondering: If I want to monitor that my velcity field is indeed divergent free, what should I check ?

I would say that in order to verify that the continuity equation is well fulfilled, I should verify that

"flow = flow_tools.GlobalFlowProperty(fchannel.solver)

flow.add_property('dx(u)+vy', name='divu_eq')",

is well close to 0.

But in order to verify that my field is well divergent free, I should verify that

"f_div_op = de.operators.differentiate(fchannel.u,'x')+de.operators.differentiate(fchannel.v,'y')"

is well close to 0.

In practice, my continuity equation is indeed verified (residual of the order 1e-17), but my field is not always divergent free as the maximum of the second operator (above) in my domain can be of the order of 1e-5 for highly transient events (turbulence for example). Is it normal, or is it an indication that my time step is too large, or should I not verify the divergence free property of my velocity field using the "second operator" ?

Thank you in advance,

Pierre-Yves

[Daniel Lecoanet]

Hi Pierre-Yves,

The divergence is dx(u) + dy(v). As you saw above, vy is not exactly equal to dy(v) — they are discrepant by the truncation error of the simulation. If the discrepancy is getting larger than is comfortable, that means that the truncation error in the simulation is too large. A large truncation error typically means you need higher spatial resolution.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6cd9eacc-42f4-49df-93ae-34089d8ecffdn%40googlegroups.com.

[Yves Dubief]

Hi Daniel,

Following on this thread, could the truncation error on  vy be responsible for the graph below? The graph shows RMS u', v', w' of velocity fluctuations in a 3D turbulent channel flow simulated with d2. The wall-normal direction is z. The asymptotic behaviors as z goes to zero are u' ~ z, v' ~ z and w' ~ z^2. The first two components exhibit the expected behavior in this compensated graph (u'_i/z) but w' does not. Although it does not appear to modify statistics in the bulk, this issue can't really be overlooked. 

Incompressibility is enforced using 

self.problem.add_equation('dx(u) + dy(v) + wz = 0')

self.problem.add_equation('wz - dz(w) = 0')

which follows d2 examples. 

Do you have any suggestion how to fix this problem? I am guessing, but I haven't tried, that d3 solves this problem with the tau formulation of the divergence of the velocity field?

Thanks

Yves

[Pierre-Yves Goffin]

Hi everyone,

By the way, I was also wondering: why is it needed to enforce the continuity equation in the form "dx(u) + dy(v) + wz = 0" and why not like "dx(u) + dy(v) + dz(w) = 0" ? I just tried the second way of doing and I'm getting errors.

Thank you in advance
Pierre-Yves

[Yves Dubief]

I am also surprised it does not work. I had tried also. I’ll ask Miguel.

Yves

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/ukpie5xcNzs/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/23206714-136b-450b-be34-f8465ea6fb8cn%40googlegroups.com.

[Daniel Lecoanet]

Hi,

Any time there is a vertical derivative, d2 automatically adds in a new tau term to the equation. If you use

problem.add_equation("dx(u) + dy(v) + wz = 0”)

problem.add_equation("wz - dz(w) = 0”)

Then you’re actually solving the equations

dx(u) + dy(v) + dz(w) + t1(x,y)*P(z) = 0

wz - dz(w) + t1(x,y)*P(z) = 0

Whereas if you use

problem.add_equation("dx(u) + dy(v) + dz(w) = 0”)

problem.add_equation("wz - dz(w) = 0”)

Then you’re actually solving the equations

dx(u) + dy(v) + dz(w) + t2(x,y)*P(z) = 0

wz - dz(w) + t1(x,y)*P(z) = 0

Which is a different set of equations, and will give different results.

In d3 we allow you to specify exactly how you want to introduce taus into your equations. If you’re worried about incompressibility, then you can enforce it exactly in d3 by not including a tau term in the incompressibility equation. Personally, I would be surprised if it makes a significant difference, but please let us know if it does!

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/F6DAD377-6AFE-4CCE-BA86-403AA429FA5C%40gmail.com.

[Pierre-Yves Goffin]

Hi,

Thank you. I would like to test to get rid of the tau term in the incompressibility equation as you mentioned. My problem is the "rayleigh_benard" 2d test case in d3. The equations and bcs are then:

problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())

problem.add_equation("trace(grad_u) + tau_p= 0")

problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")

problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - u@grad(u)")

problem.add_equation("b(z=0) = Lz")

problem.add_equation("u(z=0) = 0")

problem.add_equation("b(z=Lz) = 0")

problem.add_equation("u(z=Lz) = 0")

problem.add_equation("integ(p) = 0") # Pressure gauge

How could I "remove" the "tau_p" variable in the first equation but still ensure that the system is "well posed" ? Sorry it is quite dark the way the tau-terms work for me.

Pierre-Yves

[Calum Skene]

Hi Pierre-Yves,

This tau term is because div(constant)=0. Therefore, there is one equation in which the LHS is always zero so the problem is not solvable. The tau_p is a constant term which absorbs this redundancy. As adding tau_p adds one degree of freedom to the problem you need to add one extra equation which is that the integral of pressure is zero. This sets the gauge for pressure (as otherwise adding a constant to pressure doesn't change the equations). There's some information in the gauge conditions page here

https://dedalus-project.readthedocs.io/en/latest/pages/changes_from_v2.html#gauge-conditions

Hopefully this helps explain a bit why this tau term is there. I'm not sure if there is any benefit to trying to remove this tau_p term.

Best,

Calum

[Calum Skene]

A quick followup as I just realised that you may want to remove all tau terms from the continuity equation.

You could try

lift = lambda A, n: d3.Lift(A, lift_basis, n)

problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())

problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("dt(b) - kappa*lap(b) + lift(tau_b1,-1) + lift(tau_b2,-2) = - u@grad(b)")
problem.add_equation("dt(u) - nu*lap(u) + grad(p) - b*ez + lift(tau_u1,-1) + lift(tau_u2,-2) = - u@grad(u)")

problem.add_equation("b(z=0) = Lz")
problem.add_equation("u(z=0) = 0")
problem.add_equation("b(z=Lz) = 0")
problem.add_equation("u(z=Lz) = 0")
problem.add_equation("integ(p) = 0") # Pressure gauge

Which keeps tau_p for reasons explained before, but gets rid of the tau_u1 term from the continuity equation. This is just another way of introducing tau's to the equation (hopefully no typos).

Best,

Calum

[Pierre-Yves Goffin]

Hi Calum,

Thank you for the answer. When you say that the divergence of a constant is 0, are you speaking of the "0-mode" of my fields ? In the sense that for the first (constant) mode, in my divergence equation, I have a 0 = 0 such that I need to have another constant term in order to fix that indetermination?

I actually have no clue whether it could help or not. I'm just quite surprised, because even for the "Rayleigh-Benard" test case, the maximum divergence can go as high as 1e-2, which is quite a large error (for a full Fourier fluid simulation, the maximum divergence is typically of the order of machine precision). In fact, I have the same kind of issue for my own simulation and was trying to understand where it comes from and if it was possible to get rid of this error. Daniel then suggested to remove all the tau-terms (although he was quite not convinced that it would change anything). This is why I'm currently trying to get rid of it for the continuity equation.

Best,

Pierre-Yves

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/817fc6ab-d026-48ea-a25c-799a8c451763n%40googlegroups.com.

[Calum Skene]

Hi Pierre-Yves,

Exactly, the 0 mode of the fields. After discretisation the equation for the 0th mode becomes tau_p=0 which is solvable. Without it you get a row of zeros and you'll get a matrix is exactly singular error.

If you use the formulation I sent I think you will get machine precision for div(u)=0, but you'll maybe get more error in the continuity equation to compensate for enforcing this. I'm not 100% sure without trying. Something to maybe try is to look at how large the tau terms are in the original formulation. If they aren't small the resolution is perhaps not high enough as Daniel suggested. Increasing the resolution should make the tau terms smaller and then div(u) should also become smaller.

Best,

Calum

[Daniel Lecoanet]

Hi,

My experience is that simulations in which div(u)=0 to machine precision are no more accurate than simulations which have tau terms in the div(u) equation. If div(u) is of order 10^-2, that means that the error in everything else is also 10^-2. If you set div(u)=0 to machine precision, the error in everything else will still be 10^-2, but you won’t have any easy way to realize that. We pick parameters of the test problems so that the simulations are marginally resolved so people can run interesting-looking simulations quickly on a local machine. You can see Gibbs phenomenon in the KH example.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/dfacb610-7843-4fb0-af17-55bf1925cad0n%40googlegroups.com.

[Pierre-Yves Goffin]

Hi,

Thank you for the code suggestion Calum ! I exactly copy pasted, what you suggested, but it still doesn't work (I'm getting a "Factor is exactly singular error"). Right above what you said, I have (I don't know if there could be an issue here):

p = dist.Field(name='p', bases=(xbasis,zbasis))

b = dist.Field(name='b', bases=(xbasis,zbasis))

u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))

tau_p = dist.Field(name='tau_p')

tau_b1 = dist.Field(name='tau_b1', bases=xbasis)

tau_b2 = dist.Field(name='tau_b2', bases=xbasis)

tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)

tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)

# Substitutions

kappa = (Rayleigh * Prandtl)**(-1/2)

nu = (Rayleigh / Prandtl)**(-1/2)

x, z = dist.local_grids(xbasis, zbasis)

ex, ez = coords.unit_vector_fields(dist)

lift_basis = zbasis.derivative_basis(1)

Also, how can I have access to the tau-terms ? I thought that the tau-terms where not actually computed in Dedalus (in the sense that the matrices where manipulated in order to remove the corresponding row), but perhaps it changed in d3 ? And how can I assess whether the tau-terms are too large or not ? Should I compare it with the other coefficients or so ?

Thank you for your answer as well Daniel ! I'm not sure to perfectly grasp what you're saying though. Do you mean that since the solver solves everthing in a quite "abstract" way, the errors of each equations should be comparable, such that if I have a divergence of that order of magnitude, it simply means that I did not included enough modes in the formulation of the problem ? Then is there any general and "clean" way of assessing the accuracy of the solution (or do I have to compute the maximum residual in the grid for each equation) ? Sorry I don't exactly understand what example you are speaking about (when mentioning the "KH example"). Is it the Kelvin-Helmoltz instability and in which version of Dedalus ? Also, I though that the Gibbs phenomenon was the fact that the spectral representation of a field presents oscillations around a discontinuity, but here since the solution should be smooth, I don't quite see the issue.

Does it seems normal then that problems that only involve periodic boundary conditions typically produce more accurate solutions ? How could I understand this result ? Does it seem logical that a Chebyshev basis typically induces greater error ?

Sorry for all these questions and thank you for your help !

Pierre-Yves

[Calum Skene]

Hi Pierre-Yves,

I get the same singular error. Not sure why as of yet, I mostly do spherical shell simulations and a similar thing works there. I think it is an issue with my code rather than your part and I should have tried it before sending!

As for access to the tau terms. In Dedalus v3 you have control over how the tau terms work. In Dedalus v2 it automatically handles the tau terms in a manner similar to what the v3 Rayleigh-Bénard script does. There's a good description here https://dedalus-project.readthedocs.io/en/latest/pages/tau_method.html

In the v3 Rayleigh-Bénard script the tau terms are created as fields that only depend on the xbasis, and are solved for along with the other physical variables (u,b,p). They can be accessed in the same way that you would access u or b. Attached is a modified example that tracks the maximum tau term (this is based on code written by Ben Brown here https://github.com/bpbrown/shell_benchmark/blob/main/c2001_case1.py ) as well as maximum of the divergence over the run. I just ran this a few times for different z resolutions and if the resolution is increased the maximum of the tau terms and divergence decreases (they are fairly similar in magnitude so a high tau term gives a high divergence). Maybe this helps systematically find a resolution where everything is resolved. Hope this helps.

I am not sure about the Kelvin-Helmholtz question or the others about periodic vs Chebyshev so I'll leave these for someone else to answer :)

Best,

Calum

[Pierre-Yves Goffin]

Thank you very much for the help Calum !

[Daniel Lecoanet]

Hi Pierre-Yves,

Here is the KH example problem: https://dedalus-project.readthedocs.io/en/latest/pages/examples/ivp_2d_shear_flow.html

There is no general way of knowing the errors in each of your variables, or knowing if certain variables are more or less accurate than others. I’m just reporting what I have found. The way to determine the accuracy of your solution is to compare to a high-resolution benchmark solution. For instance, you can look here:

https://ui.adsabs.harvard.edu/abs/2016MNRAS.455.4274L

That example is compressible, but notice that the equations are conservative, so the total mass, momentum, and energy should all be conserved. This is not imposed directly in Dedalus, so we find large errors in, e.g., the total mass (~1e-5 for the Delta rho/rho = 1 case). Athena is a conservative code, so total mass, momentum and energy are conserved to machine precision. Despite the fact that Dedalus has large errors in the total mass, it gives much more accurate solutions than Athena.

So along similar lines, I’m suggesting that you can make sure div(u)=0 to machine precision, but I don’t think it will make your solution any more accurate. Don’t necessarily think it will make it less accurate though.

Hope that helps,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/36a0abbd-813b-489d-8b8a-4a97dc2d8188n%40googlegroups.com.

[Calum Skene]

Hi Pierre-Yves,

No worries, just in case it is useful, the example that is singular can be fixed in this way (also inspired by Ben Brown's code https://github.com/bpbrown/shell_benchmark/blob/main/c2001_case1.py ).

A tau term needs to be added to the divergence equation though, so it doesn't achieve what I thought! 

Best,

Calum

lift_basis = zbasis.derivative_basis(1)

lift = lambda A,n: d3.Lift(A, lift_basis, n)

# Problem

problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())

problem.add_equation("div(u) + tau_p + ez@lift(tau_u2,-1) = 0")

problem.add_equation("dt(b) - kappa*lap(b) + lift(tau_b1,-1) + lift(tau_b2,-2) = - u@grad(b)")

problem.add_equation("dt(u) - nu*lap(u) + grad(p) - b*ez + lift(tau_u1,-1) + lift(tau_u2,-2) = - u@grad(u)")

problem.add_equation("b(z=0) = Lz")

problem.add_equation("u(z=0) = 0")

problem.add_equation("b(z=Lz) = 0")

problem.add_equation("u(z=Lz) = 0")

problem.add_equation("integ(p) = 0") # Pressure gauge

[Pierre-Yves Goffin]

Thank you to both of you !

So if I'm right Daniel, in order to check the accuracy of a solution, you rather suggest to refer to the physics, than on any numerical "verification" of the equations (in the sense that there is no general way to verify the accuracy to which an equation has been solved) ? And what about the solution that Calum proposed: check the magnitude / relative magnitude of the tau-terms to verify that they don't exceed a certain threshold ?

Calum, I updated the code based on your suggestion and now it runs, but it doesn't necessarily make the divergence smaller. Thank you for the help though !

Best,

Pierre-Yves