pecos / tps

Torch Plasma Simulator
BSD 3-Clause "New" or "Revised" License
8 stars 2 forks source link

Input parsing for rxns appears broken #178

Closed trevilo closed 2 years ago

trevilo commented 2 years ago

When using multiple Arrhenius reactions, it appears that only the final reaction gets its correct parameter values. To show this, a new branch called chem-read-bug has been created where, for Arrhenius rxns, the parameters are dumped to the screen from the Chemistry class constructor (see here).

A sample input file to recreate the issue is as follows:

#---------------------
# TPS runtime controls
#---------------------

# choice of solver
[solver]
type = flow

#---------------------------------------
# Controls for flow solver
#---------------------------------------
[flow]
mesh = meshes/inline-rectangle-quad.mesh # irrelevant as long as it exists
order = 1
integrationRule = 0
basisType = 0
maxIters = 100
outputFreq = 10
useRoe = 0
enableSummationByParts = 0
fluid = user_defined
refLength = 1.
equation_system = navier-stokes
axisymmetric = True
timingFreq = 10

[io]
outdirBase = output-plasma
enableRestart = False
restartMode = singleFileReadWrite

[time]
cfl = 0.5 # p=1
integrator = rk3

[initialConditions]
rho = 1.632853e+00
rhoU = 0.0
rhoV = 1.1558
rhoW = 0.
pressure = 101325.

[boundaryConditions/wall1]
patch = 2
type = inviscid

[boundaryConditions/wall2]
patch = 3
type = viscous_isothermal
temperature = 300.

[boundaryConditions/inlet1]
patch = 4
type = subsonic
density = 1.2
uvw = '0.0 20.0 15.0'
mass_fraction/species3 = 1e-7
mass_fraction/species1 = 1e-12
mass_fraction/species2 = 0.9999999

[boundaryConditions/wall3]
patch = 5
type = viscous_isothermal
temperature = 300.

[boundaryConditions/outlet1]
patch = 6
type = subsonicPressure
pressure = 100000.

[boundaryConditions]
numWalls = 3
numInlets = 1
numOutlets = 1

#--------------------------------------------------------------
[plasma_models]
ambipolar = True
two_temperature = False
gas_model = perfect_mixture
transport_model = argon_minimal
chemistry_model = n/a
const_plasma_conductivity = 50.0

[atoms]
numAtoms = 2

[atoms/atom1]
name = 'Ar'
mass = 3.99e-2

[atoms/atom2]
name = 'E'
mass = 5.478e-7

[species]
numSpecies = 3
background_index = 2

[species/species3]
name = 'Ar.+1'
composition = '{Ar : 1, E : -1}'
formation_energy = 1.521e4
initialMassFraction = 1.00e-6
perfect_mixture/constant_molar_cv = 1.5

[species/species1]
name = 'E'
composition = '{E : 1}'
formation_energy = 0.0
initialMassFraction = 1.37e-11
# if gas_model is perfect_mixture, requires heat capacities from input.
# molar heat capacities in unit of universal gas constant.
perfect_mixture/constant_molar_cv = 1.5

[species/species2]
name = 'Ar'
composition = '{Ar : 1}'
formation_energy = 0.0
initialMassFraction = 0.9999989999863
perfect_mixture/constant_molar_cv = 1.5

[reactions]
number_of_reactions = 2

[reactions/reaction1]
equation = 'Ar + E <=> Ar.+1 + 2 E'
reaction_energy = 1.521e4
reactant_stoichiometry = '1 1 0'
product_stoichiometry = '2 0 1'
model = arrhenius
arrhenius/A = 7.279635e+08
arrhenius/b = 0.0
arrhenius/E = 1.101111e+06
detailed_balance = False

[reactions/reaction2]
equation = 'Ar.+1 + 2 E <=> Ar + 2 E'
reaction_energy = -1.521e4
reactant_stoichiometry = '2 0 1'
product_stoichiometry = '1 1 0'
model = arrhenius
arrhenius/A = 7.279635e+08
arrhenius/b = 0.0
arrhenius/E = 1.101111e+06
detailed_balance = False

Really the only important thing here is that there is more than 1 reaction. Anything else can be changed as long as the input is valid so that the Chemistry class is constructed.

When running this input using e8bda30ca53ea5ac1d4ed37edf9fdad826a04883, I get the following output to the screen:

name, input index, mixture index: E, 0, 1
name, input index, mixture index: Ar, 1, 2
name, input index, mixture index: Ar.+1, 2, 0

Inlet mass fraction of background species will not be used. 

Inlet mass fraction of electron will not be used. 
number of elements on rank 0 = 80
min elements/partition       = 80
max elements/partition       = 80
Instantiating Arrhenius rxn with A = 1.147651e-316, b = 0.000000e+00, E = 1.000000e+00
Instantiating Arrhenius rxn with A = 7.279635e+08, b = 0.000000e+00, E = 1.101111e+06

Only the second reaction is getting the correct parameter values.

Note that with 1 rxn, the parameters are correct. With 3 rxns, the first and third reactions get the correct parameters but the second are wrong. So, something is clearly going wrong with the pipeline that parses the inputs and passes them along to the Chemistry class.

trevilo commented 2 years ago

With the changes in bf615d8, the case from above gives the following output:

name, input index, mixture index: E, 0, 1
name, input index, mixture index: Ar, 1, 2
name, input index, mixture index: Ar.+1, 2, 0

Inlet mass fraction of background species will not be used. 

Inlet mass fraction of electron will not be used. 
number of elements on rank 0 = 80
min elements/partition       = 80
max elements/partition       = 80
Instantiating Arrhenius rxn with A = 7.279635e+08, b = 0.000000e+00, E = 1.101111e+06
Instantiating Arrhenius rxn with A = 7.279635e+08, b = 0.000000e+00, E = 1.101111e+06

which are the correct values.

The 1 and 3 Arrhenius rxn cases also work. To complete, need to test using a mixture of types of reactions (e.g., Arrhenius plus tabulated).

trevilo commented 2 years ago

To check with tabulated reactions in the mix, I changed the reaction block from the input above to the following:

[reactions]
number_of_reactions = 3

[reactions/reaction3]
equation = 'Ar.+1 + 2 E <=> Ar + 2 E'
reaction_energy = -1.521e4
reactant_stoichiometry = '2 0 1'
product_stoichiometry = '1 1 0'
model = arrhenius
arrhenius/A = 7.279635e+08
arrhenius/b = 0.0
arrhenius/E = 1.101111e+06
detailed_balance = False

# a fictitious reaction only to check the tabulated value.
[reactions/reaction2]
equation = 'Ar + E <=> Ar.+1 + 2 E'
reaction_energy = 1.521e4
reactant_stoichiometry = '1 1 0'
product_stoichiometry = '2 0 1'
model = tabulated
tabulated/filename = 'ref_solns/reaction/excitation.3000K.ion1e-4.h5'
tabulated/x_log = True
tabulated/f_log = True
tabulated/order = 1

detailed_balance = False

[reactions/reaction1]
equation = 'Ar + E <=> Ar + E'
reaction_energy = 0
reactant_stoichiometry = '1 1 0'
product_stoichiometry = '1 1 0'
model = arrhenius
arrhenius/A = 42.0
arrhenius/b = 3.14159
arrhenius/E = 1.414
detailed_balance = False

This leads to the following output:

name, input index, mixture index: E, 0, 1
name, input index, mixture index: Ar, 1, 2
name, input index, mixture index: Ar.+1, 2, 0

Inlet mass fraction of background species will not be used. 

Inlet mass fraction of electron will not be used. 
number of elements on rank 0 = 80
min elements/partition       = 80
max elements/partition       = 80
Instantiating Arrhenius rxn with A = 4.200000e+01, b = 3.141590e+00, E = 1.414000e+00
Instantiating Arrhenius rxn with A = 7.279635e+08, b = 0.000000e+00, E = 1.101111e+06

which has the correct Arrhenius parameters.