ReactionMechanismGenerator / ReactionMechanismSimulator.jl

The amazing Reaction Mechanism Simulator for simulating large chemical kinetic mechanisms
https://reactionmechanismgenerator.github.io/ReactionMechanismSimulator.jl
MIT License
70 stars 33 forks source link

No sensitivity for interface reactions #223

Closed ChrisBNEU closed 4 months ago

ChrisBNEU commented 1 year ago

This may be user error, or just a fundamental misunderstanding, but when doing transitory sensitivities (trapezoidal and normalized), I find that all of the sensitivities for the interface reactions are 0. I have attached an example using the Gas-Catalyst Interface notebook.

if you look at the sensitivities returned by transitorysensitivitiesfulltrapezoidal, all of the sensitivities for CH4 are either the gas phase or surface phase. there are no adsorption/desorption reactions, no matter what time you choose.

top 20 at t=1e-6s

image

top 20 at t=0.1s

image

but if I specifically look at all of the adsorption/desorption reactions for CH4, the values are all 0.

image

I thought maybe it was because I was normalizing, but I ran into another issue when I set normalized=false. it looks like the transitorysensitivitiesfulltrapezoidal function is only set up for the Simulation object, and not the SystemSimulation object, so I get the error 'type SystemSimulation has no field domain'. Gas-Catalyst Interface.ipynb.zip

mjohnson541 commented 11 months ago

It looks like we are only normalizing out the rate coefficients associated with Domains, when we should be doing this for interfaces as well here: https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl/blob/cb70323672becaf48a367172305f017809b72c55/src/TransitorySensitivities.jl#L160 and likely in other locations in that file as well.

ChrisBNEU commented 11 months ago

I wasn't sure if I should comment this on the pr or the issue, but here seemed more appropriate.

If I calculate the transitory sensitivity explicitly, and try to find any nonzero components of dSdt for parameter indices corresponding to the surface reactions, I get all 0s. I tried this at a number of times (1e-20, 1e-10, 1e-5, 1e-1), but it always returns a matrix that's all 0s.

image

I'd expect some of these values to be non-zero, especially at earlier times when it is mostly adsorption reactions taking place.

image image

Any ideas on why that might be? I am a little stuck. If I look at the parameter jacobian it is the same, unsurprisingly. Perhaps I am just fundamentally misunderstanding how the jacobian matrix is calculated.

mjohnson541 commented 11 months ago

Are those columns in just Jp (which is transitory sensitivities with tau = 0) entirely zero? Something should be at least minimally sensitive to a given rate coefficient. If entire columns are zero that might indicate something is wrong with the parameter Jacobian construction.

ChrisBNEU commented 11 months ago

if I do the same for Jp it also says that all elements are 0. At least, I assume that is what findall(!iszero, A) is giving me. maybe it rounds to 0 for some reason? I couldn't find anything in the documentation saying that though. it returns non-zero values for reactions it participates in, as I'd expect, just not the interface reactions apparently. image

mjohnson541 commented 11 months ago

So what it should be doing in a multi-domain simulation to calculate Jp is directly running automatic forward mode differentiation of the derivative function with respect to the parameters. However, the dispatch (jacobianp) we use to avoid feeding everything in pulls the Jp function off the solution object directly. You might instead try using jacobianpforwarddiff! and see if the issue is related to using that function we pull off the solution object.

ChrisBNEU commented 11 months ago

would I give it the same args as Jacobianp? because it's giving me an error MethodError: no method matching jacobianpforwarddiff!(<long nested list of data types used>) when I do that. Apologies, but I don't really know my way around julia that well, and there aren't really any examples using jacobianpforwarddiff! for a SystemSimulation. jacobianp I was able to copy code from the plotting functions.

mjohnson541 commented 11 months ago

No, the arguments are different, I believe the dispatch you want is here: https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl/blob/cb70323672becaf48a367172305f017809b72c55/src/Reactor.jl#L811

mjohnson541 commented 11 months ago

You'll need to make an empty Jacobian matrix to feed into it.

ChrisBNEU commented 11 months ago

not sure if I did this correctly, but I am getting that those columns are also all 0:

image