Zapdos chemistry is currently very simple. There are only two reaction types implemented:
\begin{equation} A \rightarrow B \\ 2A \rightarrow 2B \end{equation}For the first equation, there are two corresponding kernels, one corresponding to the sink for A and the other corresponding to the source for B. Before introducing those kernels, let's introduce a toy differential equation set that is the basis for our finite element problem:
\begin{equation} -\nabla^2A = -2k_1A^2\\ -\nabla^2B = 2k_1A^2 \end{equation}where the first equation is the balance/conservation equation for A, and the second is the balance equation for B. In this toy diffusion-reaction problem I have assumed that the only reaction occuring is $2A \rightarrow 2B$. Now to create the weak form in MOOSE we move all terms to the LHS and multiply by the test function $\psi$. In MOOSE, kernels are pieces of the governing/balance equation in weak form. So a reaction kernel should have the general form:
\begin{equation} -\psi m k C^c D^d \end{equation}where $m$ is an integer representing the stoichiometric coefficient for the variable that the kernel belongs to (positive for a product, negative for a reactant), $k$ is the reaction rate coefficient, and $C$ and $D$ are the reactants that determine the reaction rate with $c$ and $d$ their respective stoichiometric coefficients (if the reaction equation is elementary).
Now with that background let's introduce the kernels for the reaction systems $A \rightarrow B$ and $2A \rightarrow 2B$.
For $A \rightarrow B$ the corresponding kernel residual for $A$ the reactant is (in it's actual C++ form):
-_test[_i][_qp] * (-1.) * _reaction_coeff[_qp] * std::exp(_u[_qp])where _test corresponds to $\psi$, -1 is the stoichiometric coefficient for $A$ in the $A \rightarrow B$ reaction, and std::exp(_u[_qp]) represents the concentration of $A$. I should take a quick digression and say that the variable value $u$ represents the log of the actual concentration such that exp(_u[_qp]) represents the actual concentration. This logarithmic formulation is done to increase the stability of solution; it prevents concentrations from ever dropping below zero. The correspondoing kernel residual for $B$ the product is:
-_test[_i][_qp] * (1.) * _reaction_coeff[_qp] * std::exp(_v[_qp])In MOOSE, the variable that the kernel applies to (e.g. the variable for which we're adding a term to it's weak governing equation) is always denoted by $u$. It's common practice to denote a coupled variable by $v$. So in the two above residuals exp(_u[_qp]) and exp(_v[_qp]) both represent the concentration of $A$. Recognizing this, the only difference in the numerical values of the residuals comes from the switch in sign of the stoichiometric coeffient $m$. For the $A$ kernel, its value is -1; for the $B$ kernel, its value is +1.
With this background, hopefully it's fairly straightforward to understand the residual formulations for the reaction $2A \rightarrow 2B$. The residual for $A$ is:
-_test[_i][_qp] * (-2.) * _reaction_coeff[_qp] * std::exp(_u[_qp]) * std::exp(_u[_qp])The residual for $B$ is:
-_test[_i][_qp] * (2.) * _reaction_coeff[_qp] * std::exp(_v[_qp]) * std::exp(_v[_qp])Finally, hopefully we can understand how we might write the residuals for the yet unwritten kernels for the reaction $A + B \rightarrow C + D$.
For $A$ it should look like:
-_test[_i][_qp] * (-1.) * _reaction_coeff[_qp] * std::exp(_u[_qp]) * std::exp(_v[_qp])where exp(_v[_qp]) represents the concentration of $B$. The $B$ kernel residual would be the exact same except exp(_u) would represent the concentration of $B$, and exp(_v) would represent the concentration of $A$.
For $C$ and $D$, the kernel residual will look like:
-_test[_i][_qp] * (1.) * _reaction_coeff[_qp] * std::exp(_v[_qp]) * std::exp(_w[_qp])where exp(_v[_qp]) represents the concentration of $A$ and exp(_w[_qp]) represents the concentration of $B$ (recall that _u is reserved for $C$ (or $D$) which actually plays no role in determining its own rate of production).
Notation note: _qp represents the quadrature points at which the weak integral forms are converted into their discreet summation forms.
In [ ]: