search

mutationpp / Mutationpp

The MUlticomponent Thermodynamic And Transport library for IONized gases in C++
GNU Lesser General Public License v3.0
106 stars 57 forks source link

ChemNonEqTTv model setState() does not return given vibronic energies #139

Closed jbscoggi closed 4 years ago

jbscoggi commented 4 years ago

This is based on an email from @AlirezaMazaheri, who first reported this issue. I'm putting it here and listing as a bug for now.

I am trying to verify that if I set the thermodynamic state with a set of inputs, I can retrieve the inputs back after the state is set. I have verified this process for CNEQ and TCEQ and so far everything looks good. The issue that I’m having is with TCNEQ. Here, I set the TCNEQ state with species densities, rho_e and rho_eV. Then, I request the energies (total and vibrational) per unit mass and compare them with the inputs. The test fails for vibrational energy. No issue with total energy. Below are some data. Any hints as to why the returned vibrational energy is different? If it helps, I think this has to do with the presence of ionized species and/or electron in the mixture. Btw, I have checked and the Newton solver does converge.

air_5 - RRHO 
input rho_e: 446879.66295001
input rho_eV: 115166.618873251
expected e: 7710902.57631189
returned e: 7710902.57631189
expected eV: 1987198.460347526
returned eV: 1987198.460347526

air_11 - RRHO 
input rho_e: 261020770201.9375
input rho_eV: 26008381797.23112
expected e: 49265859.78793068
returned e: 49265859.78792006
expected eV: 4908901.65538193
returned eV: 4943141.489022819

CO2_8 - RRHO 
input rho_e: 18424484798.43302
input rho_eV: 55439895833.38162
expected e: 3477493.70305752
returned e: 3477493.703053357
expected eV: 10463895.7760787
returned eV: 10496633.90288009
jbscoggi commented 4 years ago

Hi @AlirezaMazaheri, can you provide some more details on how you compute the input species densities and mixture energies before the call to setState()?

AlirezaMazaheri commented 4 years ago

Hi @jbscoggi, here are the steps I take for TCNEQ. Similar steps are taken for CNEQ and TCEQ.

jbscoggi commented 4 years ago

Thanks @AlirezaMazaheri . We have a similar test already built into the test suite which seems to be functioning properly, meaning it doesn't have the same behavior. Could you provide the species densities you used that gave the results you initially reported so I can try the exact same set of numbers? I know they were random and maybe it's not possible to provide the same conditions as above. In that case, please provide all the data needed to reproduce the behavior for a case that does fail.

AlirezaMazaheri commented 4 years ago

Hi @jbscoggi can you point me to the test? When I run my test, it generates random species densities every time, and it's unlikely to get the same values as before. Currently, the Newton solver doesn't converge for mixtures other than 5-species air.

jbscoggi commented 4 years ago

Hey @AlirezaMazaheri, you beat me to it! I was just merging the updated test :). Please have a look at the following: https://github.com/mutationpp/Mutationpp/blob/a655ceb891b0371c9cc7e743cb52544ca2f00384/tests/test_set_state.cpp#L38-L96.

Note that our procedure is very similar to yours, except the densities are not obtained randomly. Instead, we equilibrate the mixture (for a grid of T,P values) and then for each T,P point, we use the equilibrium densities and a perturbed temperature vector which is the equilibrium temperature + a random perturbation in the range of +- 500K. We do this in order to have reasonable species densities for a given set of temperatures. Note that it is very easy to create non-physical densities/temperatures combos if you select them at random, independently, which could seriously impact the performance of the Newton solver.

AlirezaMazaheri commented 4 years ago

Hi @jbscoggi, Here are two sets of example (air_11 and CO2 mixtures) where the returned vibrational energies noticeably differ the input values. The species densities are in the same order as provided in the M++ mixture files. To provide a point of need access, here are the taken steps again:

air_5 - RRHO
rhos: {0.0831209815439599, 0.03252762384728387, 0.0005480571702340367, 0.02473956358534552, 8.688041095193619e-06}
input T: 8042.055877561177
input Tv: 14387.97157416015
input rho_e: 4877971.072365317
input rho_eV: 673333.7559139833
expected e: 34609060.57143456
returned e: 34609060.57143456
expected eV: 4777283.095268293
returned eV: 4777283.095268293
returned T: 8042.055877561177
returned Tv: 14387.97157416015
air_11 - RRHO 
rhos: {3.772083171225316e-07, 0.008217402454188135, 0.00158203164044893, 3.349905442259246e-06, 3.818030661535437e-06, 5.411162803271429e-08, 0.03034510868639386, 0.01013019524997852, 1.718604852728191e-06, 1.311304142186017e-05, 6.967273753567727e-08}
input T: 11742.36111897272
input Tv: 5078.056874545492
input rho_e: 2976361.926591469
input rho_eV: 80226.49707629484
expected e: 59175453.94295745
returned e: 59175453.94297452
expected eV: 1595047.746152252
returned eV: 2336397.410620507
returned T: 11076.42304153593
returned Tv: 6539.047348525341
CO2_8 - RRHO 
rhos: {1.634350644243097e-07, 0.003106515000745763, 0.0006160141919509619, 0.02351932703393298, 0.07072164609821735, 0.0007229419681532173, 3.470025039233278e-08, 1.019784785257949e-05}
input T: 15212.22216121357
input Tv: 9718.179857642906
input rho_e: 4375416.313420466
input rho_eV: 138236.0407565047
expected e: 44331878.31725346
returned e: 44331878.31725392
expected eV: 1400612.627209754
returned eV: 1672067.877769741
returned T: 14925.69634400333
returned Tv: 10843.8814584733
jbscoggi commented 4 years ago

Hey @AlirezaMazaheri, please take a look at #142 which is a test I'm adding that uses your specific inputs. While the test passes, I will note that the vibronic energy densities I get do not match the ones you have provided as "input" or "expected". These are the values I get:

air_5 rho_e: 4.87797e+06 rho_eV: 673334

air_11 rho_e: 2.97636e+06 rho_eV: 51269.8

CO2_8 rho_e: 4.37542e+06 rho_eV: 114226

Do you see anything wrong with the test code that I have added or perhaps if I am making different assumptions than you? Maybe you could try pulling this update and seeing if the test works for you.

AlirezaMazaheri commented 4 years ago

@jbscoggi Your air_5 values match the one I provided. The rho_e of air_11 and CO2_8 also match the input rho_e that I provided. Your air_11 and CO2_8 provided rho_eV are different than mine. I'll take a look at the test you just added to see if we are doing something different.

AlirezaMazaheri commented 4 years ago

Hi @jbscoggi. Mystery solved! I had a typo in my code and was incorrectly comparing hV with eV! ... and couldn't see the typo. You may close this issue. All is good.

jbscoggi commented 4 years ago

OK, good to hear @AlirezaMazaheri!