## 2 Mathematical Modelling and Numerical Methods
### Weighted Least-Squares Approximation  
For the first fractional step, different types of methods (such as finite-difference, finite-element, finite-volume, pseudo-spectral, and meshless) can be applied to discretise the spatial derivatives appearing in Equation (3). Here, for the sake of convenience, we adopted the meshless method [31], which is based on the weighted least-squares approximation. In the meshless context, the key is to approximate the spatial derivatives of a function \(f(\mathbf{x})\) in the computational domain \(\Omega\) discretised by a cloud of points (note: the points can be distributed arbitrarily, which resembles the use of unstructured grids for mesh-based body-conforming methods). Knowing discrete function values at the points \(\mathbf{x}_{i}\in\Omega\), \(i=1,\ldots,N\) (\(N\) is the total number of points), derivatives of a function at a given point \(\mathbf{x}\) can be estimated by the weighted least-squares (WLS) method, which relies on the discrete function values at the neighbour points being in the support domain of \(\mathbf{x}\). Here, the derivation is made for a two-dimensional (2D) domain. Its extension to a three-dimensional (3D) domain is straightforward. The 2D Taylor expansion of \(f(\mathbf{x}_{i})\) around \(\mathbf{x}\) is given as  
\[f(\mathbf{x}_{i})=f(\mathbf{x})+\sum_{\alpha=1}^{2}f_{\alpha}(\mathbf{x})(x_{i\alpha}-x_{ \alpha})+\frac{1}{2}\sum_{\alpha,\beta=1}^{2}f_{\alpha\beta}(\mathbf{x})(x_{i \alpha}-x_{\alpha})(x_{i\beta}-x_{\beta})+e_{i}, \tag{9}\]  
where \(e_{i}\) is the truncation error in the Taylor's series expansion, \(f_{\alpha}\) the derivative with respect to \(x_{\alpha}\) (the \(\alpha\)-th component of the position vector \(\mathbf{x}\)), and \(f_{\alpha\beta}\) the derivative with respect to \(x_{\alpha}\) and \(x_{\beta}\). The symbols \(x_{i\alpha}\) and \(x_{i\beta}\) denote the \(\alpha\)-th and \(\beta\)-th components of the position vector \(\mathbf{x}_{i}\), respectively. From the known function values \(f(\mathbf{x})\) and \(f(\mathbf{x}_{i})\) (\(i=1,2,\ldots,n\)), the unknowns \(f_{\alpha}\) and \(f_{\alpha\beta}\) for \(\alpha,\beta=1,2\) (note that \(f_{\alpha\beta}=f_{\beta\alpha}\)) are determined by minimising the error \(e_{i}\) for \(i=1,2,\ldots,n\). Here, \(n\) is the number of neighbouring points inside the support domain of \(\mathbf{x}\) (a 2D disk). Applying Equation (9) repeatedly for \(i=1,2,\ldots,n\), the system of equations for the five unknowns can be written as  
\[\mathbf{e}=Ma-\mathbf{b} \tag{10}\]  
with  
\[\mathbf{e}=[e_{1},e_{2},\ldots,e_{n}]^{\mathrm{T}},\]  
\[\mathbf{a}=[f_{1},f_{2},f_{11},f_{12},f_{22}]^{\mathrm{T}},\]  
\[\mathbf{b}=[f(\mathbf{x}_{1})-f(\mathbf{x}),f(\mathbf{x}_{2})-f(\mathbf{x}),\ldots,f(\mathbf{x}_{n})-f (\mathbf{x})]^{\mathrm{T}},\]  
\[M=[\mathbf{q}_{1},\mathbf{q}_{2},\ldots,\mathbf{q}_{n}]^{\mathrm{T}}\]  
where \(\mathbf{a}\) is the vector containing the five unknowns and \(M\) is a matrix in which the vector \(\mathbf{q}_{i}\) is defined as  
\[\mathbf{q}_{i}=\left[x_{i1}-x_{1},x_{i2}-x_{2},\frac{(x_{i1}-x_{1})^{2}}{2},(x_{i1 }-x_{1})(x_{i2}-x_{2}),\frac{(x_{i2}-x_{2})^{2}}{2}\right]^{\mathrm{T}}. \tag{11}\]  
For \(n>5\), the system (10) is over-determined. Hence, the unknowns in \(\mathbf{a}\) are determined by minimising the quadratic form  
\[J=\sum_{i=1}^{n}w_{i}e_{i}^{2}, \tag{12}\]where \(w_{i}=w(x_{i}-x)\) is the weight for the error at point \(x_{i}\). Standard minimisation of \(J\) with respect to \(a\) gives  
\[\mathbf{a}=C^{-1}A\mathbf{b}, \tag{13}\]  
where  
\[C=\sum_{i=1}^{n}w_{i}\mathbf{q}_{i}\mathbf{q}_{i}^{\mathrm{T}}, \tag{14}\]  
\[\mathbf{e}=Ma-\mathbf{b} \tag{10}\]0  
The equations above hold formally when \(f(\mathbf{x})\) is unknown and needs to be determined together with other unknowns \(f_{a}\) and \(f_{a\beta}\).