Uniform scalars don't stay uniform

[John Forbes]

Hi all,

With a lot of help from Sam and Nathan, I've found a stripped-down test case that exhibits non-uniformities in the metallicity despite a uniform initial metallicity. (Previous thread: [1]). This test case has:

- No radiative cooling

- Uniform grid ( MaximumRefinementLevel=0 )

- No gravity

- No multispecies

- No particle creation or feedback

The full parameter file is here [2]

That is, nothing fancy, just pure hydro, but I still get non-uniformities in the metallicity even after one time step.

I turned on grid::CheckDebug() and modified it to test whether |MetalDensity / Density - 1| was larger than some threshold (note that for ease of comparison I set the MetalDensity = Density initially, so that if they were evolved identically they should stay identical throughout the simulation). When this condition is no longer true, I tell the simulation to exit, which it does after calling DebugCheck from line 568 of Grid_SolveHydroEquations.C after a number of timesteps depending on the threshold I set (e.g. a threshold of 10^-4 ends the simulation during the second timestep; the first step is artificially shortened by setting Initialdt = 1.0e-6 in the paramter file). Plots of the Density [3], Metal_Fraction [4], and total velocity [5] are below.

So it looks like the Hydro solver is taking a uniform metallicity and turning it into a non-uniform metallicity. Is this behavior expected? Admittedly the initial conditions are quite violent and very artificial, so I don't expect the solution to be accurate, but I would have expected tracers to remain constant even in this situation.

Thank you for your help,

John

[1] https://groups.google.com/d/topic/enzo-dev/0Pf-QsB78Tw/discussion

[2] http://pastebin.com/XVdUNWnt 

[3] https://www.dropbox.com/s/2wj5d0yghtfna59/Density.png

[4] https://www.dropbox.com/s/v3qbm7i5hhuftk2/Metal_Fraction.png

[5] https://www.dropbox.com/s/2rdljxypw7k478a/velocity_magnitude.png

[John Forbes]

Hi everyone,

Good news! I was able to hack together some changes in the two-shock Riemann solver which solve the problem for this test case, i.e. I now get constant Metal_Fraction after several timesteps to machine precision. Basically the density and the color density were being treated slightly differently, so I was able to multiply the color field by the appropriate ratios of various interpolations / solutions of the density to get the color field to be effectively interpolated / corrected in the same way.

The big caveats are:

- I'm not sure this is the correct thing to do physically (maybe there's some drawback I haven't considered)

- I still don't know if the problem is solved in a more complicated situation (AMR, gravity, etc.) - I will check tomorrow. Turning on gravity may reveal hidden density floors like the one near line 138 of euler.F, and of course there may be separate issues with the flux correction, etc.

- These fixes will almost certainly be non-effective if certain PPM options (e.g. diffusion, steepening) are turned on. 

- Other solvers may have similar problems which these exact fixes won't solve.

Here are most or maybe all of the changes I made, i.e. a trimmed-down diff (the whole diff is a bit difficult to follow because of all the debugging stuff I added to understand what was happening).

http://pastebin.com/BQ0aSkD8

Thanks,

John

[John Forbes]

Hi folks,

A quick update for those still following the saga..

It turns out that the bug persists in AMR. My guess is that it's because, according to comments and the docs, the species fields are not flux-corrected directly, which means that the density and metal density are probably being treated differently, which will lead to this sort of problem. Right now my code in grid::DebugCheck() which checks the uniformity of the metal fraction stops the code in ProjectSolutionToParentGrid, though I think the problem is introduced to the parent grid earlier than this routine - working on narrowing it down now.

Also, I'm getting (potentially related?) warnings about flux correction, e.g.

P(0) -- CFRFr warn: 4.193201e-01 3.947483e-07 -8.723823e-01 -8.723827e-01 2 14 70 0 [1]

P(0) -- CFRFl warn: 7.789653e+01 7.189858e-06 7.812829e+01 -7.812828e+01 74 31 16 0 [0]

From the code it looks like these occur when some strictly positive quantity wants to become negative as a result of flux correction - I would have guessed that this should never occur. Should I take this as a sign that something else has gone wrong, or can reasonable simulations display this warning?

Thanks for reading,

John

[John Wise]

Hi John,

I'm following your saga! I was planning on looking into this problem,
but I haven't had the time to devise a test problem that could better
isolate it. But I can give advise on trying to troubleshoot this
problem. I would use two tests that need to be modified to follow a
color field:

1. 1D shock tube (ProblemType 1)
2. 2D Kelvin-Helmholtz Instability (ProblemType 8)

You can use the CollapseTest or something else that uses a color or
metallicity field as an example. These two tests are quick to run, and
the shock (boundary in the KH test) are well-defined and are easily
inspected in a debugger. Also, you can setup static refine regions
around the boundary layer to test whether AMR affects your solution.

About your changes to the Riemann solver: I'm not completely certain
multiplying the color left/right states by "dls(i,j)/dlX(i)" is the
thing to do, and I'd have to think about a better solution later.

About the CFRF warnings: this would indicate that a field is negative
that's not supposed be negative. Your warnings are saying that density
and total energy are going to be negative with the flux correction. Did
you see these warnings before you changes to the Riemann solver?

Good luck with this problem, and let us know if you have any other
questions.

Cheers,
John

> --
> You received this message because you are subscribed to the Google
> Groups "enzo-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to enzo-dev+u...@googlegroups.com.
> To post to this group, send email to enzo...@googlegroups.com.
> Visit this group at http://groups.google.com/group/enzo-dev.
> For more options, visit https://groups.google.com/groups/opt_out.

--
John Wise
Assistant Professor of Physics
Center for Relativistic Astrophysics, Georgia Tech
http://cosmo.gatech.edu

[John Forbes]

Hi John,

Thanks for the pointers - I'm sure they'll come in handy.

I was able to solve the problem for my test case by including both the changes to the Riemann solver and a very small change in the flux correction (FC) code. In particular, when the species fields are divided by density before FC happens, then multiplied by density after FC (so that species fraction is unaffected by FC), the metal field was, for whatever reason, not included in the list of fields considered to be species fields. Including it kept metallicity constant to ~machine precision.

I agree that these exact changes to the Riemann solver may not be the best option, but I can't quite put my finger on why. It seems like somehow you need to treat density and <color>_density's in exactly the same way, or you'll get the problem I saw. I guess the question is do these modifications treat colors correctly when they're not everywhere-constant. I think they might, since the advection equations being solved are the same, so whichever regime you find yourself in in the Riemann solver will be the same regime for density and for <color>_density, so the correction factors will be the density ratios I used, and the interpolation already done by the solver on the color fields will be correct except for these factors even when the color is not constant. I could definitely be convinced otherwise though.

I'm not sure about the CFRF warnings - I searched through the old shell output files and didn't see any, but maybe that's because these were run with CONFIG_OPT = high instead of debug (though they were also run with the -d option). 

I'm currently in the process of constructing a PR with these changes; I had to get a fresh fork of enzo so that it wouldn't include my new star formation / feedback routines and the AgoraRestart problem type which I did not write. I just have to confirm that the changes I made in the fresh version actually do solve the problem, and then maybe subsequent discussion can happen there. 

Thank you,

John

an email to enzo-dev+unsubscribe@googlegroups.com.

To post to this group, send email to enzo...@googlegroups.com.
Visit this group at http://groups.google.com/group/enzo-dev.
For more options, visit https://groups.google.com/groups/opt_out.

--
John Wise
Assistant Professor of Physics
Center for Relativistic Astrophysics, Georgia Tech
http://cosmo.gatech.edu

--
You received this message because you are subscribed to the Google Groups "enzo-dev" group.

To unsubscribe from this group and stop receiving emails from it, send an email to enzo-dev+unsubscribe@googlegroups.com.