Functions
static void	kernel_eval (double u, double &W)
	Computes the kernel function.

static void	kernel_deval (double u, double &W, double &dW_dx)
	Computes the kernel function and its derivative.

Function Documentation

◆ kernel_deval()

static void exaparty::kernels::legacy::kernelHydro::kernel_deval	(	double	u,
		double &	W,
		double &	dW_dx
	)

static

Computes the kernel function and its derivative.

The kernel function needs to be mutliplied by \(h^{-d}\) and the gradient by \(h^{-(d+1)}\), where \(d\) is the dimensionality of the problem.

Returns 0 if \(u > \gamma = H/h\).

Implementation of polynomial evaluation

We had originally the following calculation to compute w and dw_dx:

double w     = coeffs[0] * x + coeffs[1];
double dw_dx = coeffs[0];
 
for (int k = 2; k <= kernel_degree; k++) {
  dw_dx = dw_dx * x + w;
  w     = x * w + coeffs[k];
}

This variant does not vectorise for multiple reasons: We accumulate in a complicated way in w, which at least confuses the Intel 2023.2 compiler, and then we have dependencies between dw_dx and w, i.e. w feeds into dw_dx. So the idea to fuse two evaluation of Horner's method and to use the same idea of Horner (re-use of lower-order polynomials) between two quantities is contraproductive here.

We keep the plain Horner scheme for the polynomial w:

double w     = coeffs[0];
for (int k = 1 ; k <= kernel_degree; k++) {
  w = w * x + coeffs[k];
}

This part vectorises just fine and evaluates

\( w = c_{p} + x \left( c_{p-1} + x \left( \left( c_{p-2} + x \left( ... \right) \right) \right) \right) \)

For the derivative, we can start from Horner's scheme, which naturally leads to the vanilla implementation above. However, it is better to start from the plain polynomial write-up

\begin{eqnarray*} w & = & c_{0}x^p + c_{1}x^{p-1} + c_{2}x^{p-2} + ... + c_{p} \\ \frac{\partial }{ \partial x } w & = & pc_{0}x^{p-1} + (p-1)c_{1}x^{p-2} + (p-2)c_{2}x^{p-3} + ... + c_{p-1} \\ & = & (((((( p \cdot c_{0} ) x + (p-1)c_{1} ) x + (p-2)c_{2} ) x ... ) x + (p-(p-1)) c_{p-1} \end{eqnarray*}

which we then have rewritten into Horner's scheme again.

Implementation

We originally had marked this routine with InlineMethod, but that yields

kernels/legacy/kernel_hydro.h:196:68: warning: GCC does not allow
'always_inline' attribute in this position on a function definition
[-Wgcc-compat]

Parameters

u	The ratio of the distance to the smoothing length \(u = x/h\).
W	(return) The value of the kernel function \(W(x,h)\).
dW_dx	(return) The norm of the gradient of \(\|\nabla W(x,h)\|\).

Definition at line 251 of file kernel_hydro.h.

References kernel_gamma_inv, and kernel_ivals_f.