sph_kernel.py

import matplotlib.pyplot as plt
import numpy as np


class CubicKernel:
    def __init__(self, dim):
        self.dim = dim
        if dim == 1:
            self.sigma = 2. / 3.
        elif dim == 1:
            self.sigma = 10. * 0.5 * 1. / np.pi

    def get_wij(self, rij, h):
        frac_1_h = 1. / h
        q = rij * frac_1_h
        if q >= 0. and q < 1.:
            tmp = 1. - 1.5 * q**2. + 0.75 * q**3.
        elif q >= 1. and q < 2.:
            tmp = 0.25 * (2. - q)**3.
        else:
            tmp = 0.
        return tmp * self.sigma * frac_1_h**self.dim

    def get_dwij(self, xij, dwij, rij, h):
        if rij > 1e-12:
            frac_1_h = 1. / h
            q = rij * frac_1_h
            if q > 2.:
                tmp = 0.
            elif q > 1. and q <= 2.:
                tmp = - 0.75 * frac_1_h * (2. - q)**(2.)
            else:
                tmp = - 0.75 * frac_1_h * (4. * q - 3. * q**(2.))

            dw_drij = tmp * self.sigma * frac_1_h**(self.dim)

            # xij is x_i - x_j
            dwij[0] = dw_drij * xij[0] * 1. / rij
            dwij[1] = dw_drij * xij[1] * 1. / rij
        else:
            dwij[0] = 0.
            dwij[1] = 0.



def get_wij(rng, h):
    ck = CubicKernel(1)
    wij = np.zeros_like(rng)
    for i in range(len(rng)):
        rij = abs(rng[i])
        wij[i] = ck.get_wij(rij, 1.)
    return wij


def get_dwij(rng, h):
    ck = CubicKernel(1)
    dwij_array = np.zeros_like(rng)
    dwij = [0., 0.]
    for i in range(len(rng)):
        rij = abs(rng[i])
        xij = [-rng[i], 0.]
        ck.get_dwij(xij, dwij, rij, 1.)
        dwij_array[i] = dwij[0]
    return dwij_array


def sin_aproximation(pos, kernel):
    spacing = pos[2] - pos[1]
    h = np.ones_like(pos) * spacing * 1.2
    sin_actual = np.sin(pos)
    sin_sph = np.zeros_like(pos)
    for i in range(len(pos)):
        for j in range(len(pos)):
            # position between particles
            rij = abs(pos[i] - pos[j])
            sin_sph[i] += sin_actual[j] * kernel.get_wij(rij, h[i]) * spacing
    return sin_sph, sin_actual


def sin_deriv_aproximation(pos, kernel):
    spacing = pos[2] - pos[1]
    h = np.ones_like(pos) * spacing * 1.2
    sin_actual = np.sin(pos)
    sin_deriv_actual = np.cos(pos)
    sin_deriv_sph = np.zeros_like(pos)
    for i in range(len(pos)):
        dwij = [0., 0.]
        for j in range(len(pos)):
            # position between particles
            rij = abs(pos[i] - pos[j])
            xij = [(pos[i] - pos[j]), 0.]
            kernel.get_dwij(xij, dwij, rij, h[i])
            sin_deriv_sph[i] += sin_actual[j] * dwij[0] * spacing

    return sin_deriv_sph, sin_deriv_actual


def plot_kernel():
    rng = np.linspace(-3., 3., 100)
    wij = get_wij(rng, 1.)
    plt.plot(rng, wij)
    plt.show()


def plot_deriv_kernel():
    rng = np.linspace(-3., 3., 100)
    dwij = get_dwij(rng, 1.)
    plt.plot(rng, dwij)
    plt.show()

def plot_sin():
    ck = CubicKernel(1)
    x = np.linspace(-np.pi, np.pi, 100)
    sin_sph, sin_actual = sin_aproximation(x, ck)
    plt.plot(x, sin_actual)
    plt.plot(x, sin_sph)
    plt.show()


def plot_deriv_sin():
    ck = CubicKernel(1)
    x = np.linspace(-np.pi, np.pi, 1000)
    sin_deriv_sph, sin_actual = sin_deriv_aproximation(x, ck)
    plt.plot(x, sin_actual)
    plt.plot(x, sin_deriv_sph)
    plt.show()

def main():
    # plot_kernel()
    # plot_deriv_kernel()
    # plot_sin()
    plot_deriv_sin()


if __name__ == '__main__':
    main()