Functions

SPICY accepts input data of varying dimension, which are all converted into lists of 1D arrays.

spicy_vki.spicy.Spicy.__init__(self, data, points, basis='gauss', ST=None, model='laminar', verbose=2)

Initialization of an instance of the spicy class.

Parameters:

  • data (list of 1D numpy.ndarray) – Input data for the regression or a Poisson problem.

    For a regression, it has to have specific shapes

    • If model = ‘laminar’ and data = [X_G, Y_G], it is [u, v], i.e. the components of a 2D velocity field

    • If model = ‘laminar’ and data = [X_G, Y_G, Z_G], it is [u, v, w], i.e. the components of a 3D velocity field

    • If model = ‘scalar’, data can be [u1, u2, …, uN]. Multiple regressions can be done at once in this way and they share computations. However, they must share the same coordinates and the constraints have to be in the same data points.

    If the instance is to be used to solve the Poisson equation, this list contains the forcing term on the RHS of the Poisson equation. model must be ‘scalar’ in that case.

  • points (list of 1D numpy.ndarray) – Is a list of arrays containing the points: [X_G ,Y_G] in 2D and [X_G, Y_G, Z_G] in 3D.

  • basis (str, default ‘gauss’) – This defines the basis. Currently, the two options are

    • 'gauss', i.e. Gaussian RBFs exp(-c_r**2*d(x))

    • 'c4', i.e. C4 RBFs (1+d(x+)/c_r)**5(1-d(x+)/c_r)**5

  • ST (list of 1D numpy.ndarray, default None) – Is a list of arrays collecting Reynolds stresses. This is empty if the model is ‘scalar’ or ‘laminar’. If the model is RANSI (isotropic), it contains [uu’]. If the model is RANSA (anisotropic), it contains [uu, vv, uv] in 2D and [uu, vv, ww, uv, uw, vw] in 3D.

  • model (str, default ‘laminar’) – Must be one of ‘laminar’ or ‘scalar’. See input data on how to use it.

    Added in version 1.1.0.

  • verbose (int, default 2) – Sets the verbosity to print for the user

    • 0: no information printed

    • 1: updates on loops and multiple regressions

    • 2: details of each step (should only be used for a small number of regressions)

    Added in version 1.1.0.

General attributes:

  • X_G, Y_G, Z_G: coordinates of the point in which the data is available

    • if model = laminar, we assign the two (three) velocities u, v, (w) in 2D (3D) Note that each of these is of size (n_p, 1) to handle multidimensional inputs

    • if model = scalar, we assign the input data as u stacked on top of each other. I.e., if data = [u1, u2, u3], u is of size (n_p, 3)

If constraints are assigned:

  • X_D, Y_D, Z_D: coordinates of the points with Dirichlet (D) conditions

  • c_D / c_D_X, c_D_Y, c_D_Z: values of the D conditions for scalar/laminar

  • X_N, Y_N, Z_N: coordinates of the points with Neumann (N) conditions

  • n_x, n_y, n_z: normal versors where N conditions are introduced

  • c_N / c_N_X, c_N_Y, c_N_Z: values of the N conditions for scalar/laminar

  • X_Div, Y_Div, Z_Div: coordinates of the points with Div conditions

If clustering is done:

  • r_mM: vector collecting minimum (m) and maximum (M) radius of the RBFs

  • eps_l: scalar controlling the value of an RBF at the closest RBF neighbor

  • X_C, Y_C, Z_C : coordinates of the cluster centers/collocations

  • c_k: shape parameters of the RBFs

  • d_k: diameters of the rbfs

If problem is assembled:

  • A: matrix A in the linear system

  • B: matrix B in the linear system

  • b_1: vector b_1 in the linear systems

  • b_2: vector b_2 in the linear system

If computation is done:

  • weights: weights of the RBF regression

  • lambda: Lagrange multipliers of the RBF regression

spicy_vki.spicy.Spicy.collocation(self, n_K, Areas=None, bounds=None, method='clustering', r_mM=[0.01, 0.3], eps_l=0.7)

This function defines the collocation of a set of RBFs using the multi- level clustering first introduced in the article. Note that we modified the slightly original formulation to ease the programming; see video tutorials for more. The function must be run before the constraint definition.

Changed in version 1.1.0: Renamed to collocation, clustering is deprecated

Parameters:

  • n_K (list) – This contains the n_k vector in eq (33) in the paper; this is the list of expected particles per RBF at each level. For example, if n_K=[4,10], it means that the clustering will try to have a first level with RBFs whose size seeks to embrace 4 points, while the second level seeks to embrace 10 points, etc. The length of this vector automatically defines the number of levels.

  • Areas (list, default None) – List of the refinement regions for each clustering level. If no refinement is needed, then this should be a list of empty lists (default option). Currently not implemented in 3D. .. versionchanged:: 1.1.0

    Moved to keyword arguments and initialized as None.

  • bounds (list, default None) – Limits of the regression domain for semi-random clustering. In case it is not desired to use min-max rescaling, these give the limits.

    • if dimension is ‘2D’, it has to be [x_min, x_max, y_min, y_max].

    • if dimension is ‘3D’, it has to be [x_min, x_max, y_min, y_max, z_min, z_max].

    Added in version 1.1.0.

    Does not work with Areas.

  • method (str) – default=’clustering’ Which method to use for placing the collocation points. Before, only clustering was possible and this was also the original name of the function. At the moment, regular and semirandom require a rectangular domain.

    Added in version 1.1.0.

  • r_mM (list of two float values, default=[0.01, 0.3].) – This contains the minimum and the maximum RBF’s radiuses. This is defined as the distance from the collocation point at which the RBF value is 0.5.

  • eps_l (float, default=0.7.) – This is the value that a RBF will have at its closest neighbour. It is used to define the shape factor from the clustering results.

Assigns:

  • r_mM (list of length 2) – Minimum and maximum RBF radius.

  • els_l (float) – Value of an RBF at its closes neighbor.

  • X_C, Y_C, (Z_C) (1D numpy.ndarrays) – Collocation points of the RBFs in 2D (3D).

  • c_k, d_k (1D numpy.ndarrays) – Shape parameter and diamaeter of each RBF.

  • clust_list (1D numpy.ndarray) – Contains the clustering level of each RBF.

spicy_vki.spicy.Spicy.scalar_constraints(self, DIR=[], NEU=[], extra_RBF=True)

This functions sets the boundary conditions for a scalar problem. The function must be run after the clustering is carried out.

Parameters:

  • DIR (list of 1D numpy.ndarray) – This contains the info for the Dirichlet conditions:

    • If dimension = ‘2D’, it has to be [X_D, Y_D, c_D1, …, c_Dn]

    • If dimension = ‘3D’, it has to be [X_D, Y_D, Z_D, c_D1, …, c_Dn]

    Here X_D, Y_D, Z_D are the coordinates of the points where the value c_D1, …, c_Dn are set. This allows to reuse constraints for a scalar regression, assuming they are set in the same points.

    Changed in version 1.1.0: Since multiple regressions can now share computations for scalar regressions, the input can now accept a constraint value for every individual regression. The API remains the same in the case of regressing a single scalar

  • NEU (list of 1D numpy.ndarray) – This contains the info for the Neumann conditions:

    • If dimension = ‘2D’, it has to be [X_N, Y_N, n_x, n_y, c_N]

    • If dimension = ‘3D’, it has to be [X_N, Y_N, Z_N, n_x, n_y, n_z, c_N]

    Here X_N, Y_N, Z_N are the coordinates of the points where the values cN1, …, c_Nn are set for the directional derivative along the normal direction n_x, n_y, n_z.

    Changed in version 1.1.0: Since multiple regressions can now share computations for scalar regressions, the input can now accept a constraint value for every individual regression. The API remains the same in the case of regressing a single scalar

  • extra_RBF (bool, default=True.) – This is a flag to put extra collocation points where a constraint is set. It can improve the solution of the linear system as constraints remove degrees of freedom.

Assigns:

  • X_D, Y_D, (Z_D) (1D numpy.ndarrays) – Dirichlet constraints in 2D (3D).

  • X_N, Y_N, (Z_N) (1D numpy.ndarrays) – Neumann constraints in 2D (3D).

  • c_D (2D numpy.ndarray of size (n_D, n_q)) – Dirichlet constraints to be applied for the n_q quantities.

  • c_N (2D numpy.ndarray of size (n_N, n_q)) – Neumann constraints to be applied for the n_q quantities.

  • n_x, n_y, (n_z) (1D numpy.ndarrays) – Normals of the Neumann constraints to be applied.

spicy_vki.spicy.Spicy.vector_constraints(self, DIR=[], NEU=[], DIV=[], extra_RBF=True)

# This functions sets the boundary conditions for a laminar problem. The function must be run after the clustering was carried out.

Parameters:

  • DIR (list of 1D numpy.ndarray, default=[].) – This contains the info for the Dirichlet conditions. There are two options:

    • If dimension = ‘2D’, it has to be [X_D, Y_D, c_D_X, c_D_Y]

    • If dimension = ‘3D’, it has to be [X_D, Y_D, Z_D, c_D_X, c_D_Y, c_D_Z]

    Here X_D, Y_D, Z_D are the coordinates of the points where the values c_D_X, c_D_Y, c_D_Z are set.

  • NEU (list of 1D numpy.ndarray, default=[].) – This contains the info for the Dirichlet conditions. There are two options:

    • If dimension = ‘2D’, it has to be [X_N, Y_N, n_x, n_y, c_N_X, c_N_Y]

    • If dimension = ‘3D’, it has to be [X_N, Y_N, Z_N, n_x, n_y, n_z, c_N_X, c_N_Y, c_N_Z]

    Here X_N, Y_N, Z_N are the coordinates of the points where the values c_N_X, c_N_Y, c_N_Z are set for the directional derivative along the normal direction n_x, n_y, n_z.

  • DIV (list of 1D numpy.ndarray, default=[].) – This contains the info for the divergence-free conditions. There are two options:

    • If dimension = ‘2D’, it has to be [X_Div, Y_Div]

    • If dimension = ‘3D’, it has to be [X_Div, Y_Div, Z_Div]

    Here X_Div, Y_Div, Z_Div are the coordinates of the points where the divergence-free condition is imposed.

  • extra_RBF (bool, default=True.) – This is a flag to put extra collocation points where a constraint is set. It can improve the solution of the linear system as constraints remove degrees of freedom

Assigns:

  • X_D, Y_D, (Z_D) (1D numpy.ndarrays) – Dirichlet constraints in 2D (3D).

  • X_N, Y_N, (Z_N) (1D numpy.ndarrays) – Neumann constraints in 2D (3D).

  • c_D_X, c_D_Y, (c_D_Z) (dd numpy.ndarrays) – Dirichlet constraints to be applied for the three velocities.

  • c_N_X, c_N_Y, (c_N_Z) (dd numpy.ndarrays) – Dirichlet constraints to be applied for the three velocities.

  • n_x, n_y, (n_z) (1D numpy.ndarrays) – Normals of the Neumann constraints to be applied.

spicy_vki.spicy.Spicy.plot_RBFs(self, level=0, plot_skip=1)

Utility function to check the spreading of the RBFs after the clustering. This function generates several plots. It produces no new variable in SPICY.

Parameters:

  • level (int, default 0) – This defines the cluster level of RBF that will be visualized.

  • plot_skip (int, default 1) – Only every plot_skip’th point will be plotted to save cost for large datasets.

    Added in version 1.1.0.

spicy_vki.spicy.Spicy.Assembly_Poisson(self, n_hb=0)

This function assembly the matrices A, B, b_1, b_2 for the Poisson problem. These are eqs. (31a) - (31d) in the original paper (see also video tutorial 1 for more info)

Parameters:

n_hb (int, default 0) – When solving the Poisson equation, global basis elements such as polynomials or series expansions can be of great help. This is evident if one note that the eigenfunctions of the Laplace operator are harmonics. In a non-homogeneous problem, once could homogenize the basis. This will be proposed for the next relase (which will align with Manuel’s paper). The idea is the following: if the homogeneization is well done and the basis is well chosen, then we do not need constraints for these extra terms of the basis.

For the moment, we let the user introduce the number of extra_basis. These will be sine and cosine bases, which are orthogonal in [-1,1]. In 1D, they are defined as : sines_n=np.sin(2*np.pi*(n)*x); cos_n=np.cos(np.pi/2*(2*n+1)*x) Given n_hb, we will have that the first n_hb are sines the last n_hb will be cosines. This defines the basis phi_h_n, with n an index from 0 to n_hb**4 in 2D.

In 2D, assuming separation of variables, we will take phi_h_nm=phi_n(x)*phi_m(y). Similarly, in 3D will be phi_nmk=phi_n(x)*phi_m(y)*phi_k(z). For stability purposes, the largest tolerated value at the moment is 10!.

For an homogeneous problem, the chosen basis needs no constraints.

Warning: This feature is currently under development.

spicy_vki.spicy.Spicy.Assembly_Regression(self, n_hb=0, alpha_div=None)

This function assembly the matrices A, B, C, D from the paper (see video tutorial 1).

Parameters:

  • n_hb (int, default 0) – Also for a regression, the harmonic basis can improve the regression as they can model global trends which are similar to a low order polynomial. Furthermore, for homogenous problem, they automatically fulfill the boundary conditions.

    See the same entry in the function ‘Assembly_Poisson’

  • alpha_div (float, default None) – This enables a divergence free penalty in the entire flow field. The higher this parameter, the more SPICY penalizes errors in the divergence-free condition. This is particularly important to obtain good derivatives for the pressure computation.

spicy_vki.spicy.Spicy.Solve(self, K_cond=1000000000000.0)

This function solves the constrained quadratic problem A, B, b_1, b_2. The method is universal for 2D/3D problems as well as laminar/Poisson problems.

The input parameters are the class itself and the desired condition number of A which is fixed based on its largest and smallest eigenvalue

The function assigns the weights ‘w’ and the Lagrange multipliers Lambda to the class. The weights are computed for the min/max scaled problem, i.e. the right hand-side of the linear system is normalized. The assigned weights are rescaled by self.scale_U to get the real, physical quantities

Parameters:

K_cond (int, default 1e12.) – This is the regularization parameter. It fixes the condition number (see Video 1) The estimation is based such that the regularize matrix has the condition number k_cond. For this, we compute the max and the min eigenvalue.

spicy_vki.spicy.Spicy.get_sol(self, points, order=0, shape=None)

This function evaluates the solution of the linear system on an arbitrary set of points. The function is used for both the function itself as well as the analytical derivatives.

Parameters:

  • points (1D numpy.ndarray) – Contains the points at which the source term is evaluated

    • if dimension is ‘2D’, it has to be [X_P, Y_P].

    • if dimension is ‘3D’, it has to be [X_P, Y_P, Z_P].

  • order (1D numpy.ndarray, default 0) – Order of the derivative solution

    • order = 0 : Function without any derivative

    • order = 1 : First derivative

    • order = 2 : Second derivative

    • …

    Added in version 1.1.0.

  • shape (tuple, default None) – Shape in which to reshape the solution arrays. This is a QOL feature. The product of the tuple entries has to have the same length as the point array.

Returns:

U_sol – The solution always depends on model = ‘scalar’/’laminar’ and dimension = ‘2D’/’3D’ If order = 0, give the original function

  • If scalar, the solution is [U1_sol, … Un_sol].

  • If laminar and 2D, the solution is [U_sol, V_sol]

  • If laminar and 3D, the solution is [U_sol, V_sol, W_sol]

If the scalar solution is only [U1_sol], it is returned as U1_sol so that a list of length 1 does not need to be unpacked by the user.

If order == 1, give the first derivative of all quantities

  • If ‘scalar’ and ‘2D’, the output is [dU1DX, dU1DY, …, dUndx, dUndy]

  • If ‘scalar’ and ‘3D’, the output is [dU1DX, dU1DY, dU1DZ, …, dUndx, dUndy, dUndz]

  • If ‘laminar’ and ‘2D’, the solution is [dUdx, dUdY, dVdX, dVdY]

  • If ‘laminar’ and ‘3D’, the solution is [dUdx, dUdY, dUdX, dVdX, dVdY, dVdX, dWdX, dWdY, dWdX]

Return type:

list of numpy.ndarray or numpy.ndarray

spicy_vki.spicy.Spicy.get_source_term(self, points, rho)

This function evaluates the source term on the right hand side of equation (21) in the paper (see video tutorial 1 for more info).

Parameters:

  • points (list of 1D numpy.ndarrays) – Contains the points at which the source term is evaluated. If the model is 2D, then this has [X_P, Y_P]. If the model is 3D, then this has [X_P, Y_P, Z_P].

  • rho (float) – Density of the fluid. Has to be a float number.

Returns:

source_term – Normal pressure in equation (29) evaluated in the given points.

Return type:

1D numpy.ndarray

spicy_vki.spicy.Spicy.get_pressure_neumann(self, points, normals, rho, mu)

This function evaluates the Neumann boundary conditions for the pressure integration in equation (29) from Sperotto et al. (2022) (see video tutorial 1 for more info)

Parameters:

  • points (list) – Contains the points at which the Neumann constraint is evaluated. There are two options:

    • If the model is 2D, then this has [X_P, Y_P].

    • If the model is 3D, then this has [X_P, Y_P, Z_P].

  • normals (list) –

    Contains normals of the points at which the Neumann constraint is evaluated. There are two options:

    • If the model is 2D, then this has [n_x, n_y].

    • If the model is 3D, then this has [n_x, n_y, n_z].

  • rho (float) – Density of the fluid. Must be larger than 0.

  • mu (float) – Dynamic viscosity of the fluid. Must be larger than 0.

Returns:

P_neu – Normal pressure in equation (29) evaluated in the given points.

Return type:

1D numpy.ndarray