Setting a floor density

[Min-Kai Lin]

Dear colleagues,

What is the correct way to set a floor value during a simulation? At the moment, to set a minimum value of a variable Q, I simply put

while solver.proceed:

...

          Q['g'][Q['g']<floor] = floor 

...

Indeed, printing min(Q) within this loop indicates it is working. However, the hydrodynamic outputs of Q, set via "snapshots.add_task" (before the time integration loop) still display negative values. 

Thanks,

Min-Kai

[Daniel Lecoanet]

Hi Min-Kai,

After you modify Q on the grid, it gets transformed to spectral coefficients for timestepping. Then for it to be output, it is transformed back to grid space. If your floor is low, I wonder if this grid -> coefficients -> grid round-trip transform is leading to a slightly different answer which pushes you below zero again. You could test this by running Q.require_coeff_space() and Q.require_grid_space() before printing min(Q). Note that we drop the Nyquist mode, and changing only a few values of Q will likely project onto the Nyquist mode.

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/f9aed5b4-2322-4691-93bb-9976c2b62e66n%40googlegroups.com.

[Min-Kai Lin]

Dear Daniel,

Many thanks for the swift clarification. Just to confirm: this issue only affects the output values, not the simulation itself, i.e. the floor values are correctly set/used before each timestep?

Cheers,

Min-Kai

[Daniel Lecoanet]

Hi Min-Kai,

Dedalus evolves the spectral coefficients of the problem variables. If you make a change to the grid data that does not project onto the spectral modes we evolve (e.g., if you add a multiple of the Nyquist mode), then you have not actually changed the simulation state at all. If you modify the simulation data on the grid, then you need to do a grid -> coeff -> grid transform to see how you have actually modified the spectral coefficients of the problem, which is exactly what happens when you do an output. So the output are giving a correct depiction of the simulation.

I’m not sure if this is possible in your simulations, but in the past when I have solved equations for the density (which should be >0), I instead wrote the equation in terms of log(rho), which automatically imposes rho>0.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/db9b0529-9d27-471b-b33e-5f61c380e30an%40googlegroups.com.

[Min-Kai Lin]

Dear Daniel, 

Thanks for the detailed explanation. I realize the following question is about methods rather than Dedalus, which may not be appropriate here, but perhaps you/other users have some quick tips on maintaining positivity in a spectral code. 

After some experimentation, it seems that modifying the grid data cannot remove negative densities (e.g. the filter in this paper). I suppose this is because abruptly resetting selected grid values to zero doesn't play well with the spectral transformation. 

I also tried transforming the variable as you suggested, but I found that the total mass isn't conserved. More specifically, the relevant equation in my problem looks like 

d/dt(Q) + v*grad(Q) = - div(Q*F), 

where div(v)=0, Q is a density that should be >0, and F is a vector function of other problem variables. 

I tried solving for

Y = Q0*(1 + ln(Q/Q0)),

where Q0 is the initial value of Q. So the equation for Y is

d/dt(Y) + v*grad(Y) = -grad(Y)@F - Q0*div(F),

I then add a Laplacian diffusion term D*lap(Q) to the RHS. However, this formulation does not conserve the total mass

V = volume_integral(Q) = volume_integral(Q0*exp(Y/Q0 - 1))

, where I used Dedalus' volume_integral function. Non-conservation manifests in the nonlinear regime. 

Thanks very much,

Min-Kai

[Daniel Lecoanet]

Hi Min-Kai,

If you add D*lap(Q) to the RHS of your equation for Y, then you’re right that the equation does not conserve Q. If you start with the equation

d/dt(Q) + v*grad(Q) = - div(Q*F) + D*lap(Q),

and divide by Q, you get

d/dt(log(Q)) + v*grad(log Q) = - grad(log(Q)).F - div(F) + D*lap(Q)/Q

one can check

Q^{-1} lap(Q) = lap(log(Q)) + grad(log(Q)).grad(log(Q))

So defining Y = logQ, we get

d/dt(Y) + v.grad(Y) - grad(Y).F - div(F) + D*lap(Y) + D*|grad(Y)|^2

That equation should conserve the volume integral of Q.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6a87cfbd-6c3e-433d-aad6-54a70460b3e3n%40googlegroups.com.

[Min-Kai Lin]

Hi Daniel,

Many thanks for the tip! Just to follow up in case anyone else stumbles upon this thread. I think it works now after transforming the diffusive term properly. I note that fairly high resolution is needed to conserve mass; e.g. with 128 modes/box length, the error is about 0.2%, which in my case is already practical. 

Cheers,

Min-Kai