2017-10-18 06:29 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002578OpenFOAMContributionpublic2017-07-16 21:03
Reportertniemi 
Assigned Tohenry 
PrioritynormalSeverityminorReproducibilityalways
StatusclosedResolutionsuspended 
Product Versiondev 
Target VersionFixed in Version 
Summary0002578: Comments regarding chemistryModel jacobian
DescriptionNote: 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.
TagsNo tags attached.
Attached Files

-Relationships
+Relationships

-Notes

~0008213

tniemi (reporter)

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.

~0008219

henry (manager)

> 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?

~0008225

tniemi (reporter)

> 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.

~0008226

henry (manager)

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.

~0008380

henry (manager)

Have you completed the testing of this proposed change? Should I close this report for now and you can re-open when the testing is complete?

~0008401

tniemi (reporter)

No, I haven't yet done any more work for this. You can close the report for now.

~0008402

henry (manager)

Pending completion of the testing for the contribution.
+Notes

-Issue History
Date Modified Username Field Change
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
2017-07-11 09:52 henry Note Added: 0008380
2017-07-16 19:47 tniemi Note Added: 0008401
2017-07-16 21:03 henry Assigned To => henry
2017-07-16 21:03 henry Status new => closed
2017-07-16 21:03 henry Resolution open => suspended
2017-07-16 21:03 henry Note Added: 0008402
+Issue History