quokka-astro / quokka

Two-moment AMR radiation hydrodynamics (with self-gravity, particles, and chemistry) on CPUs/GPUs for astrophysics
https://quokka-astro.github.io/quokka/
MIT License
46 stars 12 forks source link

Odd-even instability error in diffusion limit #541

Closed chongchonghe closed 8 months ago

chongchonghe commented 9 months ago

Describe the bug The advecting radiation pulse test has what seems to be an odd-even instability issue when the pulse is not advecting (v = 0) and it goes away when v > 0.

Method 1:

epsilon = {S_corr * S_corr, S_corr, S_corr, S_corr} for even indices and epsilon = {S_corr, 1.0, 1.0, 1.0} for odd indices.

radhydro_pulse_velocity original

Method 2:

epsilon = {S_corr, 1.0, 1.0, 1.0} for even indices and epsilon = {1.0, 1.0, 1.0, 1.0} for odd indices.

radhydro_pulse_velocity

Method 2 seems to work well for the MarshakAsymptotic problem as well: compare method 1 and method 2:

CleanShot 2024-02-20 at 13 18 26

To Reproduce Steps to reproduce the behavior:

  1. Compile this problem 'test_radhydro_pulse'
  2. Run this problem 'test_radhydro_pulse'
  3. See the error in the velocity curve.
BenWibking commented 9 months ago

Can you fill out the sections requested in the description for this issue?

BenWibking commented 9 months ago

My proposed solution is to modify the reconstruction so that Equation 27 of [1] is obtained in the diffusion limit. This will require interpolating between the standard reconstruction and Equation 27 as a function of $\tauj$ and $\tau{j+1}$.

[1] https://ui.adsabs.harvard.edu/abs/2001JQSRT..69..475L/abstract

chongchonghe commented 9 months ago

@BenWibking @markkrumholz See PR #544 . What Ben suggested is also used in the THC_M1 code (Radice2022). So, I followed their approach and also tried a few variations. However, none of them worked. I'll have to other methods, like switching the order of flux and energy update.

BenWibking commented 9 months ago

No, I don't think it is the case that my suggestion above is the same method as that used by THC_M1. (It is a little unclear since the notation in that paper is not fully defined in the text.)

My suggestion above was to modify the reconstruction, not the Riemann solver. As far as I can tell, this is not what is described in the THC_M1 paper, nor what you have implemented in PR #544.

chongchonghe commented 8 months ago

OK. Here are the results for the advecting pulse test in the static/dynamic diffusion limit with low (64 grids)/high (512 grids) resolution. I'm using the Method 2 mentioned above for odd-even correction.

Static diffusion, low-res: static-diffusion-64grids-oddeven-type1

Static diffusion, high-res: pulse-test-static-diffusion

Dynamic diffusion, low-res: combined

Dynamic diffusion, high-res: pulse-test-dynamic-diffusion

BenWibking commented 8 months ago

Does it make any difference if you set hydro.artificial_viscosity_coefficient=0.1?

chongchonghe commented 8 months ago

With higher viscosity, the oscillation is slightly suppressed but not very much.

Advecting pulse in dynamic diffusion limit:

CleanShot 2024-03-05 at 12 14 00

CleanShot 2024-03-05 at 12 14 25

BenWibking commented 8 months ago

Reasonable values are between 0 and 1. The default is zero. The value recommended in the literature is usually 0.1.

It looks like a value of 1 suppresses a lot of the oscillations, though not completely. Does setting it to 1 at high resolution also suppress the oscillations?

BenWibking commented 8 months ago

My guess is that the oscillations and the temperature differences are not related. I think the temperature difference is caused by the flattening or "clipping" of extrema that occurs with PPM limiting (this also happens with lower-order reconstruction).

I would try to see if disabling the PPM limiters for this problem and seeing if it fixes the problem (i.e., only using Equation 22 of https://arxiv.org/pdf/2110.01792.pdf, and turning off the subsequent steps). It's possible that this is only stable when artificial viscosity is nonzero.

chongchonghe commented 8 months ago

My guess is that the oscillations and the temperature differences are not related.

I tend to believe in this as well. With 2 times better resolution (nx=128), the oscillation is at the same level, but the temperature difference is gone. See figure below:

CleanShot 2024-03-05 at 12 51 23

However, this is not better than viscosity = 0.

CleanShot 2024-03-05 at 12 58 33

I would try to see if disabling the PPM limiters for this problem and seeing if it fixes the problem (i.e., only using Equation 22 of https://arxiv.org/pdf/2110.01792.pdf, and turning off the subsequent steps). It's possible that this is only stable when artificial viscosity is nonzero.

I'll try that. I also believe the reconstruction is a big issue here because when setting the reconstruction order to 1 or 2, the solution to the static un-advecting pulse is totally wrong, as I mentioned in the this PR.

chongchonghe commented 8 months ago

At this point, I would call it a victory no matter what, since the dynamic diffusion test 'passed' with only 128 grids without artificial viscosity.

BenWibking commented 8 months ago

Yes, agreed.

I think completely fixing the oscillations would require re-doing the asymptotic analysis with the odd-even fix, and then doing a linear perturbation analysis to see what the growing modes are. Since the oscillations in velocity are very small relative to the sound speed, however, I don't think it's worth it.

chongchonghe commented 8 months ago

I would try to see if disabling the PPM limiters for this problem and seeing if it fixes the problem (i.e., only using Equation 22 of https://arxiv.org/pdf/2110.01792.pdf, and turning off the subsequent steps). It's possible that this is only stable when artificial viscosity is nonzero.

No luck with this. I tried with nx=128 and artificial_viscosity_coefficient=0.1. The magnitude of the oscillation is 3x bigger.

BenWibking commented 8 months ago

I would try to see if disabling the PPM limiters for this problem and seeing if it fixes the problem (i.e., only using Equation 22 of https://arxiv.org/pdf/2110.01792.pdf, and turning off the subsequent steps). It's possible that this is only stable when artificial viscosity is nonzero.

No luck with this. I tried with nx=128 and artificial_viscosity_coefficient=0.1. The magnitude of the oscillation is 3x bigger.

Yes, I expect it would make the oscillations worse. Instead, the question is whether it makes the low-resolution temperature agree better.