Initial guess for temperature(s) can be problematic when setting state with conserved variables

[mgoodson-cvd]

When using set_state with vars = 0, which uses conserved variables (partial densities rhoY and total internal energy density rhoe and vibrational energy rhoev for 2T), a Newton iteration is required to determine the temperature. This is done for 1T using getTfromRhoE and for 2T using solveEnergies. The Newton iteration requires an initial guess for the temperature. In both instances, the initial guess for the temperature comes from m_T, which is the stored temperature in memory. If using set_state for the first time since initializing a mixture object, the initial guess will be 300 K because this is what m_T is initialized to in the StateModel constructor. However, on all subsequent calls using the same mixture object, m_T will simply be whatever was in memory from the last time the temperature was set.

This isn't necessarily wrong or even a bad idea. Imagine using M++ in a code and looping over cells, using a single mixture object, passing in conserved vars. The first cell performs the Newton solve, gets its temperature, carries on. If the next cell in your loop isn't in too different of a state, which is usually the case, the previous solution is probably a great initial guess. However, I see three potential issues with this. The first is when the next cell isn't in a similar state, such as when crossing a shock. However, Newton does pretty well with a sane initial guess, so this is likely okay. The second is when the previous cell didn't converge for some reason. If this happens, you're probably in trouble anyway, but your next call to M++ might fail due to the bad initial guess, even if you gave it sane inputs. This seems bad. The third major issue is repeatability. The fact that my initial guess depends on the previous call means that if I change the cell order, I can get different results. As long as my Newton iteration converges, this is likely not a major issue.

Of these three issues, the second issue is the most concerning to me. There are several potential remedies.

The first is to always start from the same initial guess. There is actually commented out code in the 2T solveEnergies function that sets both m_T and m_Tv to 500 K for the initial guess, every time. Is this safer than the above? Yes. Is it less performant? Likely yes.

The second is to allow the user to specify an initial guess. This can actually already be accomplished by just calling set_state with vars=1, passing in temperatures and partial densities. The temperatures will be used to set m_T and m_Tv and as long as you don't do anything else, those will be preserved in memory and used as the initial guesses when you call set_state with vars=0. This method does not require code modification, but instead places the burden on the user to call set_state twice. Alternately, an optional input could be added to set_state with the temperature(s) array.

The third is to check whether or not the stored value for m_T and m_Tv are sane before using it as the initial guess. If the value is "bad" (which would need to be defined...), then reset to some sane value, like 500 K. I would define a "good" temperature to be between 50 K and 20,000 K. The 1T getTfromRhoE already "clamps" the temperature at 50 K, so that seems like a good lower bound. I think I've seen an upper bound on temperature somewhere in M++ but I can't recall where; but if you're using M++ above 20,000 K, you may want to think again.

I am leaning toward the option 3, at minimum. I am leaning away from option 1 because it is less performant. The advantage of option 2 is that a guess provided by the user is likely to be the best. Doesn't mean we can't still have option 3 in there.

[jbscoggi]

Hey @mgoodson-cvd, I'm actually surprised we are still using the previous temperature as the initial condition. I think I have probably commented that out and used a fixed temperature of 500K several times in the past. In fact, using the previous temperature led to a very serious and difficult-to-diagnose bug when computing finite-difference jacobians in a CFD solver. Basically, when the solver would perturb each conserved quantity to compute the finite-difference, the starting point for the temperature would be different each time leading to slightly different results for each perturbation. This in turn would lead to non-zero terms in the jacobian which should have been strictly zero. In generally, it's definitely safer to have a fixed initial guess, or at least a deterministic algorithm for choosing a good guess. For example, you could always compute the enthalpy at two fixed temperatures and do a linear interpolation with the given enthalpy to find an initial temperature. That should always give the same results when providing the same enthalpy and might improve the Newton convergence. This might be a case of premature optimization however. It is probably worth getting some statistics of number of Newton iterations required for a bunch of representative cases to see how much the initial guess matters for performance.

[mgoodson-cvd]

@jbscoggi Very interesting. Commented-out lines left in source, such as the line setting the temperature to 500 K here, are rarely without a purpose or story -- good to learn that for this one!

I can do some testing on convergence. As I mentioned in the original post, I am still somewhat hesitant to always start at a fixed temperature of say 500 K since it will hurt performance on sane cases. However, it will also solve some issues so really the best answer is to allow an initial guess, and if none is provided, then fall back to 500 K. I can't see a scenario in which using what was in memory is desired.