Avoid machine precision subtraction issue in simple transport #607
Conversation
| Ddiag[j] += Xloc[i] * dbintemp; | ||
| } | ||
| } | ||
| for (int i = 0; i < NUM_SPECIES; ++i) { |
There was a problem hiding this comment.
I think we could merge this for loop into the previous one, to avoid doing an extra N^2 loop, but also remove the if statement, i.e.
- we add term1=0 at the start of the loop
- in the first nested loop (j < i) we add term1 += yloc[j]
- Add another loop j=i+1, j < NUM_SPECIES which just adds the remaining part of term1
- Construct Ddiag in the i loop.
|
I had issues with pure species even before we added this change. I've been running pure hydrogen jets, and the diffusion solver was failing. I gradually increased the inflow mass fraction, and it was fine up to 99.9% hydrogen, but it would fail as soon as it hit 100%. Running with a fixed Lewis number ran fine, so I thought it was related to the mixture averaged formulation for pure species. |
|
@ThomasHowarth Good call on also making the changes for SRK. Unless I'm missing something your proposed changes to the loop structure would result in the same number of FLOPS and I don't think we can move the full Ddiag within the first i loop so I did not implement that. One potential issue with the trace stuff for pure component diffusion is that it permanently modifies the input Yloc after the transport function is called, which is unexpected behavior. I In general the whole process here when applied to pure components makes me a bit uneasy because we're essentially dividing by zero with a numerical hack that makes it sort of work. Not sure what a better approach would be though. |
|
I'm going to merge this so we can get PelePhysics updates in PeleC going again. @ThomasHowarth we can revisit your suggestion on rearranging the loops later if you'd like |
|
Sorry, somehow I completely missed that PR. I'll take a look. |
In simple transport coefficient, the transport coefficient for species i is roughly proportional to$\sum Y_{j \neq i}/\sum X_{j \neq i}$ . This leads to zero divided by zero numerical issues for Yi = Xi= 1. These are sidestepped by adding an epsilon 1e-15 quantity to all mass fractions just for purposes of the diffusion coefficient calculation to ensure no mass fraction is ever exactly unity in these calculations.
This leads to numerical precision issues when we compute$\sum Y_{j \neq i}$ as $1 - Y_i$ . Summing a bunch of quantities of O(1e-15) retains full floating point precision for sum, and therefore full floating point precision for the resulting diffusion coefficient. However, taking the difference between two O(1) quantities results in an O(1e-15) quantity with only O(1e-1) precision, and therefore results in a diffusion coefficient that can have O(1e-1) error, as seen in the EB-InflowBC test here: AMReX-Combustion/PeleC#971.
The proposed fix is to revert to computing$\sum Y_{j \neq i}$ instead of computing $1 - Y_i$ , which unwinds part of the performance benefit from #595.