Formulations-and-Numerical-Schemes/Implementation.md

+The aim of this section is to provide information about algorithm used in the library to compute some of the equations presented in section **Spatial discretization**.
+
+To numerically compute the NSW equations we need to evaluate several operator like physical flux, numerical flux and integrals. Here we will present features to better understand the method use to calculate those terms.
+
+# Example of cells integrals
+
+Practically, to evaluate the term $`\int_{\Omega_{e}} \widetilde{F}(W_{h}(x), z_{h}(x)) \cdot \nabla \varphi(x) dx`$ of the NSW equations wo do as follow :
+
+- Evaluate $`W_{h}`$ and $`z_{h}`$ at quadrature points on the reference element :
+
+```math
+W_{h}(\widehat{x_{i}}) = \sum_{j=1}^{p} W_{j} \widehat{\varphi_{j}}(\widehat{x_{i}})
+```
+- Evaluate the flux term $`\widetilde{F}(W_{h}, z_{h})`$ at quadrature points.
+	
+- Compute the integral using quadrature formulas :
+
+```math
+	\begin{aligned}
+		\int_{\Omega_{e}} \widetilde{F}(W_{h}(x), z_{h}(x)) \cdot \nabla \varphi(x) dx &= \int_{\widehat{\Omega_{e}}} \widetilde{F}(W_{h}(\widehat{x}), z_{h}(\widehat{x})) \cdot J_{F_{e}}^{-1} \widehat{\nabla} \widehat{\varphi}(\widehat{x}) det(J_{F_{e}}) d\widehat{x} \\
+		&= \sum_{i=1}^{n} \omega(\widehat{x_{i}}) \widetilde{F}(W_{h}(\widehat{x_{i}}), z_{h}(\widehat{x_{i}})) \cdot J_{F_{e}}^{-1} \widehat{\nabla} \widehat{\varphi}(\widehat{x}) det(J_{F_{e}}) d\widehat{x} 
+	\end{aligned}
+```
+
+Remark : as the position of the quadrature points do not change over time, we store the values of the $`\widehat{\varphi}(\widehat{x_{i}})`$ and $`\widehat{\nabla} \widehat{\varphi}(\widehat{x_{i}})`$ for efficiency.
+
+# Data storage in **Aerosol**
+
+The unknown in the library are stored following a pattern described here :
+
+- **String Vector** : Is a class containing information about data, like the name of the vector, the size of the system associate with the data, the tag in Pampa, etc... It do not contain any data or pointer.
+
+- **Scheme** : $`\texttt{map<StringVector, PampaVector>}`$
+	In spatial scheme (and in mesh) we create map that link *String Vector* to *Pampa Vector*.
+
+- **Pampa Vector** : It is basically a $`\texttt{map<int, DataHandler>}`$ that, for each type of element in the mesh, associate a *Data handler*.
+
+- **Data handler** : It is the class that actually provide the pointer of data.
+
+In order to manage memory in parallel, pointer are allocated by calling **Pampa** when we fill the $`\texttt{map<StringVector, PampaVector>}`$ in the spatial scheme. No global memory allocation is performed outside spatial scheme (except in the mesh class, but you should not have to deal with it!). 
+
+As a consequence, in **Aerosol** and so in **Uhaina**, the data is store following this pattern :
+
+- Type of element.
+- Entities in type.
+- Degree of freedom (DOF).
+- Variables of the system.
+
+*[element_type_1][element_type_2] ... [element_type_n]*
+
+$`\hspace{1cm}`$ *|*
+
+$`\hspace{0.5cm}`$ *[element_number_in_type_1] ... [element_number_in_type_n]*
+
+$`\hspace{2cm}`$ *|*
+
+$`\hspace{1.5cm}`$ *[DOF_1] ... [DOF_N]*
+
+$`\hspace{3cm}`$ *|*
+
+$`\hspace{2.5cm}`$ *[variable_1] ... [variable_n]*
+
+
+
+
+
+
+
+# Algorithm
+
+This section describe the algorithm used to compute the residual. It is located in the function :
+
+
+```c++
+	AdvectionBoussinesqOptim::fillResidualVectorT;
+```
+
+To evaluate the NSW equations we use the following algorithm :
+
+1. Preprocessing, compute interpolation and projection matrix.
+
+2. Compute residual at each time step :
+
+	(a) `Start the communications.` 
+
+
+	(b) Loop over local elements of the mesh :
+
+		- Perform interpolation at quadrature points.
+
+		- Evaluate flux and source term on element.
+
+		- Evaluate elements integrals.
+
+		- Store face quadrature points data.
+
+	(c) Loop over boundary face of the mesh and compute boundary flux. Store.
+
+	(d) Loop over internal, non-frontier face of the mesh and compute numerical flux. Store.
+
+	(e) Loop over internal elements (with no frontier faces) and add the contribution of faces integrals.
+
+	(f) `Wait for communication to finish.`
+
+	(g) Loop over overlap elements and perform interpolation at quadrature points. Store.
+
+	(h) Loop over internal, frontier faces of the mesh. So face with an element on the overlap. Compute numerical flux and store.
+
+	(i) Loop over border elements (with frontier faces) and add the contribution of faces integrals.
+