|View Issue Details|
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0002578||OpenFOAM||Contribution||public||2017-06-09 13:05||2017-06-16 17:12|
|Target Version||Fixed in Version|
|Summary||0002578: Comments regarding chemistryModel jacobian|
|Description||Note: this is not strictly a bug report, but more of a commentary/sharing of experiences|
Lately, I have been testing the EDC and TDAC with the GRI-mechanism. Generally, I have had good success in some test cases, but in others the chemistry solution tends to crash quite often.
I recently saw a presentation, where it was shown that much better stability and performance could be obtained by replacing the jacobian in the chemistryModel with an analytical version obtained from pyJac. (There were also other changes such as replacing LUSolver with Lapack) Based on this, I looked at the jacobian-function to see what might cause stability/performance problems.
If I have understood correctly, currently the jacobian does the following:
- First calculates dcdt
- Then calculates (semianalytically) d(dcdt)dc terms
- Finally calculates d(dcdt)dT numerically
However, I can see several potential issues:
1. According to ODESystem, the jacobian is supposed to return dfdx, which I think would mean d(dcdt)dt. However, here the jacobian calculates dcdt, which would correspond to dydt. (Also, dcdt contains space for dTdt which is simply ignored here.)
2. When calculating d(dcdt)/dc terms, it is assumed that the R.kf and R.kr do not have derivatives with respect to concentrations. This is not true with thirdBodyReactions or FallOffReactions.
3. d(dTdt)/dc terms and d(dTdt)/dT are assumed to be zero, which is not generally true.
I guess the jacobian doesn't need to be very accurate, so are these issues intentional optimizations or bugs?
I made some tests where I included the derivatives of R.kf and R.kr in the jacobian calculated with simple forward finite difference. I have attached a comparison of the chemFoam gri-tutorial case. As can be seen, the amount of integration steps clearly drops with high tolerances. However, the jacobian evaluation is now more costly and depending on tolerances, the overall performance may increase or decrease.
I have attached the modified chemistryModel, which I used in the tests. It is hard code optimized for GRI, but it could be easily made more general if the reaction rates would hold a flag telling whether or not they have a non-zero derivative wrt concentrations. Even better would be, if the reaction rates could calculate their own derivatives, but this would require more work.
There might also be some other optimizations possible. I tested caching the Kc value when calculating R.kr derivatives and that already shaved off 10 % of the computation time.
I also tried to address concerns 1. and 3., but they didn't seem to affect the results much.
|Tags||No tags attached.|
Related to optimizations, I have also tried to look on how to optimize the LUDecompose method. The current implementation looked rather complicated, so I tried to create a more simplified version to see if it would be faster.
I have attached alternative versions of LUDecompose and LUBacksubstitute and an artificial benhmark testcase. In the testcase the alternative implementation is roughly 30-40 % faster than the current implementation for medium to large matrices (n>10). However, the performance depends largely on compiler optimizations, so the results may vary. I only used a single compiler (gcc 4.8.3), but several different machines. On an older machine (Intel X5650), the performance of the stock LUDecompose seemed to drastically decrease when n>500 (probably cpu cache limitation) and for these large matrices the performance difference was 300-400 %. On newer machines the performance difference stayed at 30-40 %.
I also benchmarked the alternative implementation against Eigen, which is supposedly a very optimized linear algebra library. For small matrices (n<50) there were no diffences, but with n~1000 Eigen was about 3 times faster than the alternative LUDecompose. So for large matrices there is still room to the state of the art. Of course, for such a large systems using sparse matrices would be more reasonable.
I have also attached a comparison of all the chemFoam tutorials with the stock settings and with the alternative LUSolver. For h2 and gri there is practically no difference, but on larger cases the difference is noticeable.
Note that the new implementation uses different decomposition algo, so the decomposition is not identical with the stock one.
> replacing the jacobian in the chemistryModel with an analytical version
This would be a good idea and would not require interfacing to a Python library to achieve. Someone would need to spend a bit of time working out the form of the temperature derivative which would be a bit tedious but not fundamentally difficult.
> I guess the jacobian doesn't need to be very accurate, so are these issues intentional optimizations or bugs?
Simplifications rather than bugs, and not explicit optimizations. If you would like to contribute a more complete form of the evaluation of the Jacobian I would be happy to include it.
> I have attached alternative versions of LUDecompose and LUBacksubstitute
I would be happy to include these developments if you are sure that they do not adversly affect the operation of other parts of the code which rely on them. What tests have you run so far?
> Simplifications rather than bugs, and not explicit optimizations. If you would like to contribute a more complete form of the evaluation of the Jacobian I would be happy to include it.
Unfortunately, I don't currently have enough time to do this properly. Besides chemistryModel, similar changes should also be made to TDACChemistryModel, which is slightly more complex due to the mechanism reduction.
> I would be happy to include these developments if you are sure that they do not adversly affect the operation of other parts of the code which rely on them. What tests have you run so far?
I have run the chemFoam tutorials and the Matrix-Test program. I have also checked that in the benchmark tests both the new and old decompositions give similar solutions to random linear systems. However, I did spot a slight mistake in the submitted version and I have uploaded a new one. In the previous version I was not taking the abs of the diagonal element when checking for the largest pivot. This led to unneccessary pivoting and reduced stability in cases where the diagonal was negative.
After the fix I have not observed any performance degradations or other issues. However, apart from the chemFoam tests and benchmarking tests I have not done other testing. Also, I have only tried gcc compilers, so I don't know whats the performance difference with eg. clang.
Below is a list of known differences between the implementations:
- The methods produce a slightly different decompostions and pivots. Both are valid, but they do not produce exactly the same results.
- Both methods use partial pivoting so their stability should be roughly the same.
- If the stock implementation encounters a singular matrix during pivoting, it will replace the zero with SMALL. This allows the code to always continue, although the solution of the linear system will likely be garbage. The alternative implementation will halt in these situations.
|While it is likely that the new LUDecompose and LUBacksubstitute will work fine for other the purposes they are used for in OpenFOAM this will need to be tested before commiting to avoid surprises for the other users of OpenFOAM.|
|2017-06-09 13:05||tniemi||New Issue|
|2017-06-09 13:05||tniemi||File Added: comparison.png|
|2017-06-09 13:06||tniemi||File Added: chemistryModel.C|
|2017-06-13 14:21||tniemi||Note Added: 0008213|
|2017-06-13 14:21||tniemi||File Added: LUcomparison.png|
|2017-06-13 14:22||tniemi||File Added: benchmark.zip|
|2017-06-13 14:22||tniemi||File Added: alt_LUDecompose.zip|
|2017-06-14 10:50||henry||Note Added: 0008219|
|2017-06-16 13:29||tniemi||File Added: alt_LUDecompose_fixed.zip|
|2017-06-16 13:39||tniemi||Note Added: 0008225|
|2017-06-16 17:12||henry||Note Added: 0008226|