*2.2. Computational Implementation*

Next, we describe some of the topics involved in the computational implementation of the ML-enhanced, FEM-based method utilized to tackle the drop deformation problems proposed in Section 3. All code was written in the C programming language and compiled using the open-source GNU-C compiler gcc. Efficiency and robustness are two features pursued in this implementation; accordingly, all simulations presented in this work can be run in a commodity personal computer with a 4-core processor. All 2D plots were made using the open-source library matplotlib [55] and/or the Ti*k*Z and PGF packages [56], while the 3D graphics make use of Mayavi [57]. For data manipulation and data analysis, we use the pandas [58] package.

### 2.2.1. Machine Learning Enhancement

Radial basis functions are a family of kernel methods [59] widely used in Machine Learning (ML) for their flexibility in generating spaces of trial functions with excellent approximation properties. This close relation to approximation theory allows for a solid mathematical background, which proves extremely useful when offering bounds in supervised learning problems [43] for the necessary number of training data required to provide a model trained with such data and having a small generalization error. Combining these concepts with novel applications in the field of image processing [60,61], we leverage RBFs to reconstruct the discrete solution of the polymer stress tensor over the whole domain, using the scattered data provided by the values of the polymer stress tensor defined at each of the position of the ensembles. The main advantage of this idea is the mesh independence property guaranteed by the RBF approach; the relevance of this feature cannot be overstated in the context of "micro-macro" methods for non-Newtonian flows [29,34,41]. Traditional methods use particles for the kinetic modeling of the polymer solution; however, spatial mesh refinement hinders the accuracy of the stochastic model (given the reduced number of particles surrounding a certain mesh node), thus requiring compromising on the accuracy of either the "micro-" (internal configurations of the polymer and "extra-"stress tensor) or the "macro-" (velocity, pressure) scales. This is no longer the case with our approach: now, using the smooth approximation built via RBFs, we can evaluate the tensor at any position (e.g., a certain mesh point), virtually decoupling the "micro-" and "macro-"scales and paving the way for adaptive mesh refinement in "micro-macro" methods.

Thus, at each instant of time, we retrieve the positions {**<sup>x</sup>***i*} , *i* = 1, .., *Nens* of all the ensembles scattered over the domain. Then, we compute *<sup>τ</sup>p*,*i*, the value of the polymer stress tensor at the ensemble *i* using Kramers' expression averaging over the *Nd* dumbbells contained at the ensembles:

$$\mathbf{r}\_{p,i} = nk\_B \boldsymbol{\Theta} \left[ \frac{1}{N\_d} \sum\_{j=1}^{N\_d} \mathbf{F} \left( \mathbf{Q}\_j^i \right) \otimes \mathbf{Q}\_j^i - \mathbf{I} \right]. \tag{5}$$

As kernels of approximation, we focus on compactly-supported RBFs [59,62] for their computational efficiency and remarkable accuracy, since the matrices resulting from a suitable formulation of the approximation problem are sparse. In particular, we use Wendland's *<sup>ϕ</sup>*3,*<sup>k</sup>*(*C*2*<sup>k</sup>*), *k* = 0, .., 3, suitable for problems of space dimension *d* ≤ 3, and degree of smoothness 2*k*:

$$\begin{cases} \wp\_{3,0}(r) = (1-r)\_{+' \prime}^2 \\ \wp\_{3,1}(r) = (1-r)\_{+}^4 (4r+1)\_{\prime} \\ \wp\_{3,2}(r) = (1-r)\_{+}^6 (35r^2 + 18r + 3)\_{\prime} \\ \wp\_{3,3}(r) = (1-r)\_{+}^8 (32r^3 + 25r^2 + 8r + 1)\_{\prime} \end{cases} \tag{6}$$

so that outside the support size *χ*, *<sup>ϕ</sup>*3,*k*(*<sup>r</sup>* > *χ*) = 0, whereas inside the support, they are positive definite on R*<sup>d</sup>* and provide optimal convergence rate <sup>O</sup>(*h*Ω)<sup>2</sup>*k*+*d*+1, with *h*Ω the fill distance [62]. Using these kernels, we build an interpolant *s* of the form:

$$s(\mathbf{x}) = p(\mathbf{x}) + \sum\_{i} \lambda\_i q(\|\mathbf{x} - \mathbf{x}\_i\|),\tag{7}$$

with *p*(**x**) a quadratic function defined by coefficients *ci* and and *λi* a set of real-valued interpolation coefficients for the CSRBFs. The formulation can then be reduced to a system of equations equivalent to a saddle-point problem of the form [Q P; P<sup>T</sup> 0]] [l; c] = [f; 0], which can be solved efficiently using the two-step approach described in [52] to obtain the solution vectors [l;c] required to build the interpolant *s*. After building the interpolant and having reconstructed the polymer stress tensor, we can straightforwardly evaluate each component at any given point (e.g., at one mesh node).

The use of CSRBFs as approximation kernels calls for the determination of the support size *χ*, for each of the ensembles scattered over the domain. The technique used to perform this potentially expensive computation is taken from the classification problems usually found in ML: the Nearest Neighbor (NN) classification algorithm [42,43]. Using this approach, along with an implementation of kd-tree as data structures, we are able to obtain, at each instant of time, the nearest neighbors of each ensemble, within a certain support size *χ*.
