Equation based reaction constants = 0

[youareit125]

Hi,

I am working on modeling a steady-state argon plasma. Below is a section of the input file.

[ChemicalReactions] [ScalarNetwork] species = 'Ar4s Ar4p Arp Ar2p' aux_species = 'Ar e' equation_constants = 'Tgas' equation_values = '300' equation_variables = 'Te' interpolation_type = 'linear'

reactions = 'e + Ar -> e + e + Arp         : {2.34e-14*Te^(0.59)*exp(-17.44/Te)}
                    e + e + Arp -> e + Ar         : 1e-40
                    Ar2p + Ar -> Arp + Ar + Ar    : {(6.06e-6/Tgas)*exp(-15130.0/Tgas)}'

[] []

At t = 0, output prints Te = 3 and Tgas = 300 but reaction_constant0 = 0, reaction_constant1 = 1e-40 and reaction_constant2 = 0. Is there a way for crane to accept Te and Tgas into the reaction rates? Within executioner, I am using type = Steady and solve_type = 'Newton'

[smpeyres]

Double checked this against the two reaction argon and argon microdischarge problems; indeed, k = 0 for anything equation-based for t = 0. This seems to occur in all solves, but since we've almost only done Transient problems, the short dt values seems to "screen out" this effect. (Assuming these dynamical systems are sufficiently non-chaotic.) However, this completely paralyzes practical application of steady solves, since the entire system solves on t = 0 and gives steady-state output on t = 1.

@dcurreli Seems that the reaction parser has quite a few problems...

@youareit125 This is may take a bit to fully figure out, I will keep you updated.

[youareit125]

After switching the solve_type within the Executioner block from NEWTON to LINEAR, the equation-based rate constants are solved, and the calculation converges. Will try a simpler model to confirm everything makes sense in a few days.

[dcurreli]

@youareit125 Can you provide the full input file, and where/how you define Te, Tgas, etc.

Also, the Steady option is new and under-tested. I suggest to use Transient for simple problems like this. Computational time won't be an issue.

[youareit125]

@dcurreli Below is the input file. Having 'linear' for solve_type in Executioner block works whereas 'newton' won't. I would like to use Transient, but I couldn't find a way to calculate time dependent electron temperature and electron density.

https://github.com/youareit125/Crane-Steady-State/blob/06c743b576b9f7097cd0fa9ff80490f4d78231e7/Steady_TC_v1.i

[smpeyres]

@youareit125, @dcurreli and I have reviewed this issue in depth and have placed it under high-priority as fixing the initialization of equation-based reaction coefficients is crucial for both transient and steady-state problems. My hypothesis, without digging into the code, is that the reaction rate coefficients that are equation-based are updated at the end of every timestep rather than at the beginning. This should be the other way around to ensure proper initialization at t = 0.

The fix itself is expected to be quite easy (somewhere there's probably a "switch" that lets us say TIMESTEP_BEGIN instead of TIMESTEP_END or something similar), but every test across CRANE and ZAPDOS that uses the ChemicalReactions block will need to be updated to reflect this change, which will take longer. We hope to have a fix within two weeks, but we will update you if there are any delays.

[dcurreli]

Check line 370 here:

https://github.com/lcpp-org/crane/blob/devel/src/actions/AddReactions.C#370

[smpeyres]

Update: I replaced the "TIMESTEP_END" with "TIMESTEP_BEGIN" within AddReactions.C and AddZapdosReactions.C and compiled successfully, but I did not see any changes in the initial evaluation of equation-based rate coefficients. This was confirmed by the fact that all of the tests still passsed, indicating that nothing in the output changed.

Might be something deeper...

[dcurreli]

Checking for consistency between the [ChemicalReaction] block and [AuxVariables], e.g. on the Lorentz Attractor problem, replacing for example {p} with an Auxiliary Variable:

https://github.com/lcpp-org/crane/blob/devel/problems/lorentz_attractor.i

[smpeyres]

Update @youareit125 @dcurreli : I have went through all of the .C files and appended “INITIAL” to wherever “TIMESTEP_BEGIN” was set as the execution point for the evaluation of rate coefficients. This is simply a caveat of the MOOSE system that differentiates t = 0 and the beginning of all other timesteps. So far, all I have not seen any k = 0 evaluations at t = 0 with this change, but all tests still pass. This could be a nuance of the EXODIFF-type test, but not sure.

@youareit125 The rate coefficients at t = 0 had non-zero evaluations for your input file under this change, and coverged under a linear solve. I tried a netwon solved, but it failed. My assumption is that it failed due to the other numerical parameters or a lack of constraints on some variables.

[dcurreli]

@smpeyres This seems to fix the problem. Please commit your changes on a separate branch so we can test further.

[dcurreli]

Check on

• [x] Two-reaction Argon
• [x] Argon micro-discharge

[smpeyres]

PR #118 fixes this issue. Currently in devel branch.

Thank you @youareit125 for pointing this out! Let’s get in touch to discuss the steady-state numerics.