Fuego w/ FORD rates generating incorrect prefactor conversion term

[whitmanscu]

With the Fuego-generated mechanisms for kinetics with FORD rates, species production/consumption rates appear unphysically low. I believe this is due to the generated rates in comp_qfqr, where the prefactor conversion appears to not be properly implemented for FORD rates. There is conversion that happens between cgs and mks before productionRate is called, and the prefactor term seems to be intended to convert back the other way based on the mechanism units. However, the dim that is used in calling self._prefactorUnits in CPickler.py for the prefactor term is based off of the stoichiometric coefficients rather than the FORD rates, which leads to an improper conversion. If I manually adjust the prefactor terms based off of the FORD rates, my kinetics appear to give realistic production/consumption behavior.

A quick fix seems like it would be to change https://github.com/AMReX-Combustion/PelePhysics/blob/283a0475e2487ed570eb47c3118d0fc981e40e05/Support/Fuego/Pythia/pythia-0.4/packages/fuego/fuego/serialization/c/CPickler.py#L1751 to something like

if (len(reaction.ford) > 0) :                
   dim = self._phaseSpaceUnits(reaction.ford)            
else:                
   dim = self._phaseSpaceUnits(reaction.reactants)

similar to what's done in the initialization. However, I want to check and make sure this is the most appropriate fix, or if a different route might be better.

[EnnaDelfen]

Sounds like this is what is going on, and the fix proposed would solve the issue

[drummerdoc]

@EnnaDelfen Is there any additional code that would impact the analytic Jacobian that would also need fixing?

[EnnaDelfen]

I'm looking through it because he's specifically taking about productionRate_GPU I'm looking through to see if other functions are concerned

[EnnaDelfen]

Yeah so I'm thinking there are a bunch of other places where the suggested modifications need to be implemented. It's tricky to catch them all without a mechanism/test case to test the implementations. Basically wherever a dim needs to be computed or a call to _QSSsortedPhaseSpace is needed, FORD need to be taken into account. I think no one has actually ued a FORD reaction scheme in where AJ or other fancy chemistry treatment was required and that's why it simply hasn't been implemented. SHould be a matter of grepping the code for calls to _QSSsortedPhaseSpace and dim and checking results. I'm happy to review the proposed changes !

[drummerdoc]

@whitmanscu Can you work with me to create a test case that we can use to help debug the fixes? For this, I'm guessing that we can modify PelePhysics/Testing/Exec/ReactEval_C to integrate one of the FORD-containing models using both implicit and explicit integrators and that would uncover most of the implementation issues. Which model where you using, and what could we set for a set of initial conditions for the integration and a delta_t integration period?

[whitmanscu]

@drummerdoc Happy to. I'm using a 2-step global propane model from Ghani et al. (2015). I can push that model to a testing branch, or we could use something already there like the CH4 2-step model. I'm not as familiar with the expected evolution of the ReactEval_C case, but it looks like a variant of an ignition delay problem. At atmospheric pressure and phi=0.65, the adiabatic propane flame temperature is right around 1,800K and the ignition delay time for the global mechanism at, say, 1,100K is estimated at 9.5e-5s by Cantera; at 1,800K that drops to 9.2e-6s. In the PeleC case I've been running with the propane mechanism and the rk64 reactor, the time steps are on the order of 1e-7. Hopefully this gives an idea of the best way to set up the ReactEval_C case. Let me know if more information would be helpful.

[whitmanscu]

I added the C3H8 model and modified the CPickler.py comp_qfqr generation as suggested. The prefactor does appear to generate correctly for comp_qfqr after this change: https://github.com/whitmanscu/PelePhysics/tree/whitmans/FORD_prefactor_fix

Regarding ReactEval_C, I'm not set up to run on GPUs at the moment, and it looks like the more complex routine testing for things like cvode w/ analytical jacobians have that as a requirement. Let me know if there's further information or testing I can help provide.

[malihass]

closed by #414