I am looking at a tutorial create for the FENICs finite element method package. The tutorial shows a system of advection-diffusion-reaction equations occurring in a solution that is moving according to the incompressible navier-stokes equations. In the tutorial, the authors simply present the weak form, without explaining the derivation. I was hoping someone might be able to explain how to derive the weak form, or even identify some good references to look at to find these derivations.
The system of equations is below. [NOTE: I included the equations as images because I could not get the latex formatting to work for an aligned environment. I can put in the latex, but might need some direction on how to format the equations.]
The provided weak form is as below. I understand the steps in computing a weak form, especially the integration by parts piece. However, I am not clear on all of the steps in the derivation below, such as why a vector PDE is represented as a single equation. I imagine that the different species are tracked in a vector, and hence the sum makes sense. BUT, in the equation below, it is a bit difficult to see those nuanced elements. Hence I was hoping someone could explain this derivation or provide some links to resources that provide a more detailed derivation of this weak form.
This is a bit subtle, but fundamentally it isn't too difficult once you write out the spaces an inner products and stuff. For simplicity, assume that each of the species and corresponding test functions $u_i$, $v_i$ are elements of ther same function space $V$ that embeds continuously into $L^2(\Omega)$ so that we can use the $L^2$ inner product without any issues. Denote this inner product by $(~,~)$.
Consider just equation for the first species. We perform the backwards Euler time discretization to obtain the strong form $$\frac{u_1^{n+1} - u_1^n}{\Delta t} + w\cdot\nabla u_1^{n+1} - \nabla(\epsilon \nabla u_1^{n+1}) - f_1 + Ku_1^{n+1}u_2^{n+1} = 0.$$ This isn't completely decoupled from $u_2$, but bear with me for a bit. You can consider $u_2$ to be known for the time being. Since $u_1^{n+1}$ is unknown, we define $u_1^{n+1}$ as the solution of this equation; i.e., we write $\mathcal{F}_1(u_1^{n+1}) = 0$, where $$\mathcal{F}_1(\varphi):=\frac{\varphi - u_1^n}{\Delta t} + w\cdot\nabla \varphi - \nabla(\epsilon \nabla \varphi) - f_1 + K\varphi u_2^{n+1}.$$ To determine the weak form for this nonlinear equation, we multiply $\mathcal{F}_1$ by a test function from $V$ and integrate by parts, but, more precisely, we take the inner product of $F_1$ with a test function and integration by parts allows us to more easily compute a symbolic form for that inner product. Indeed, since we are using the $L^2$ inner product, the weak form is given by $$F_1(u,v) = (v,\mathcal{F}_1(u)) = \int_\Omega(\Delta t^{-1}(u - u_1^n)v + q\cdot\nabla uv + \epsilon\nabla u\cdot\nabla v - f_1v + Kuu_2^{n+1}v)~\mathrm{d}x.$$ This is standard and is likely most of what you have seen previously.
Notice that $F_1(u,v)$ takes values in $\mathbb{R}$. This should be the case with any weak form, and must also be the case when we consider all species together, so we will need to add things together. However, the way to derive that comes naturally if you pay attention to the choice of function space. Now consider all three species and their respective test functions, which all come from $V\subset L^2(\Omega)$. To consider them in a coupled manner where they all simultaneously solve the same problem, we now define a new vector function space $\mathbf{V} = V\times V\times V$ with inner product $$\langle \mathbf{u},\mathbf{v}\rangle = (u_1,v_1) + (u_2,v_2) + (u_3,v_3),$$ where $\mathbf{u} = (u_1, u_2,u_3)$ and $\mathbf{v} = (v_1,v_2,v_3)$. Similarly, after the temporal discretization we have a strong form given by $$\boldsymbol{\mathcal{F}}(\mathbf{u}) = \begin{pmatrix}\mathcal{F}_1(\mathbf{u}) \\ \mathcal{F}_2(\mathbf{u}) \\ \mathcal{F}_3(\mathbf{u}) \end{pmatrix}.$$ My definitions are slightly different than in the first example because now each PDE, denoted $\mathcal{F}_i$, depends explictly on the unknown value of each species, $u_i^{n+1}$ instead of considering one to be fixed.
We now proceed as before to derive the weak form, except now, multiplying by a test function and integrating by parts must be abstracted to instead taking the inner product of $\boldsymbol{\mathcal{F}}(\mathbf{u})$ with a test function from the vector function space $\mathbf{V}$. In particular, we have $$\mathbf{F}(\mathbf{u},\mathbf{v}) = \langle \mathbf{v}, \boldsymbol{\mathcal{F}}(\mathbf{u})\rangle = (v_1,\mathcal{F}_1(\mathbf{u})) + (v_2,\mathcal{F}_2(\mathbf{u})) + (v_3,\mathcal{F}_3(\mathbf{u})).$$ From this, we can see that the total weak form is just the sum of the weak forms of the individual PDEs with respect to different components of the vector test function $\mathbf{v}$.