# do not remove; needed to identify dependencies
from numpy import shape
import struphy.kernel_arguments.pusher_args_kernels as pusher_args_kernels
import struphy.pic.sorting_kernels as sorting_kernels
import struphy.pic.sph_smoothing_kernels as sph_smoothing_kernels
from struphy.kernel_arguments.pusher_args_kernels import MarkerArguments
[docs]
def distance(x: "float", y: "float", periodic: "bool") -> float:
r"""Return the signed one-dimensional distance ``x - y``, adjusted for periodicity on ``[0, 1]``.
Parameters
----------
x, y : float
Two coordinates on the domain.
periodic : bool
If ``True``, the result is folded into :math:`(-\tfrac{1}{2}, \tfrac{1}{2}]` so that
the shortest path across the periodic boundary is returned.
Returns
-------
float
Signed distance ``x - y``, adjusted for periodicity.
"""
d = x - y
if periodic:
if d > 0.5:
while d > 0.5:
d -= 1.0
elif d < -0.5:
while d < -0.5:
d += 1.0
return d
########################
# single-point kernels #
########################
[docs]
def naive_evaluation_kernel(
args_markers: "MarkerArguments",
eta1: "float",
eta2: "float",
eta3: "float",
holes: "bool[:]",
periodic1: "bool",
periodic2: "bool",
periodic3: "bool",
index: "int",
kernel_type: "int",
h1: "float",
h2: "float",
h3: "float",
) -> float:
r"""Perform a single-point SPH evaluation of a function :math:`\rho: [0, 1]^3 \to \mathbb R` in the following sense:
.. math::
\rho(\boldsymbol \eta_i) = \sum_{j=0}^{N-1} \rho_j\, W_h(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.
The coefficients :math:`\rho_j` must be available in the marker array, stored at some index ``self.markers[index, j]``.
In case that `derivative=k` where `k` is not zero, the `k`-th component of the gradient of :math:`\rho` is computed:
.. math::
\textrm{derivative}=k:\qquad [\nabla \rho(\boldsymbol \eta_i)]_k = \sum_{j=0}^{N-1} \rho_j \frac{\partial W_h}{\partial \eta_k}(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.
The possible choices for :math:`W_h` are listed in :ref:`smoothing_kernels`
and in :meth:`~struphy.pic.base.Particles.ker_dct`.
ATTENTION: The sum is done over all particles in the markers array (ignoring holes and ghost particles), no neighbour search is performed.
Hence, the cost of this evaluation is :math:`\mathcal{O}(N)` in the number of particles, and it should only be used for testing and verification purposes.
Boundaries are always treated as periodic (mirror images held by ghost particles are not included).
Parameters
----------
args_markers : MarkerArguments
Container holding the markers array and the total number of particles ``Np``.
eta1, eta2, eta3 : float
Evaluation point in logical space.
holes : bool[:]
1D array of length ``markers.shape[1]``. ``True`` if particle ``i`` is a hole (inactive).
periodic1, periodic2, periodic3 : bool
``True`` if the domain is periodic in that dimension.
index : int
Column index in the markers array of the coefficient :math:`\beta_k` multiplying the kernel.
kernel_type : int
Integer identifier of the smoothing kernel. See :ref:`smoothing_kernels`.
h1, h2, h3 : float
Kernel width in the respective dimension.
Returns
-------
float
SPH estimate of :math:`b` at the evaluation point.
"""
markers = args_markers.markers
n_particles = shape(markers)[1]
out = 0.0
for p in range(n_particles):
# skip holes and ghost particles (copies of markers already counted, id = -2)
if not holes[p] and markers[-1, p] != -2.0:
r1 = distance(eta1, markers[0, p], periodic1)
r2 = distance(eta2, markers[1, p], periodic2)
r3 = distance(eta3, markers[2, p], periodic3)
out += markers[index, p] * sph_smoothing_kernels.smoothing_kernel(kernel_type, r1, r2, r3, h1, h2, h3)
return out
[docs]
def box_based_kernel(
args_markers: "MarkerArguments",
eta1: "float",
eta2: "float",
eta3: "float",
loc_box: "int",
boxes: "int[:,:]",
neighbours: "int[:,:]",
holes: "bool[:]",
periodic1: "bool",
periodic2: "bool",
periodic3: "bool",
index: "int",
kernel_type: "int",
h1: "float",
h2: "float",
h3: "float",
) -> float:
r"""Perform a single-point SPH evaluation of a function :math:`\rho: [0, 1]^3 \to \mathbb R` in the following sense:
.. math::
\rho(\boldsymbol \eta_i) = \sum_{j=0}^{N-1} \rho_j\, W_h(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.
The coefficients :math:`\rho_j` must be available in the marker array, stored at some index ``self.markers[index, j]``.
In case that `derivative=k` where `k` is not zero, the `k`-th component of the gradient of :math:`\rho` is computed:
.. math::
\textrm{derivative}=k:\qquad [\nabla \rho(\boldsymbol \eta_i)]_k = \sum_{j=0}^{N-1} \rho_j \frac{\partial W_h}{\partial \eta_k}(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.
The possible choices for :math:`W_h` are listed in :ref:`smoothing_kernels`
and in :meth:`~struphy.pic.base.Particles.ker_dct`.
The sum is restricted to the 27 neighbouring boxes of the box containing
:math:`\boldsymbol\eta_i`, making the cost :math:`\mathcal{O}(1)` in the number
of particles when the kernel support is proportional to the box size.
Parameters
----------
eta1, eta2, eta3 : float
Evaluation point in logical space.
loc_box : int
Box of the evaluation point.
boxes : 2d array
Box array of the sorting boxes structure.
neighbours : 2d array
Array containing the 27 neighbouring boxes of each box.
markers : array[float]
Markers array.
Np : int
Total number of particles.
holes : bool
1D array of length markers.shape[1]. True if markers[:, i] is a hole.
periodic1, periodic2, periodic3 : bool
True if periodic in that dimension.
index : int
Column index in markers array where the value multiplying the kernel in the evaluation is stored.
kernel_type : int
Number of the smoothing kernel.
h1, h2, h3 : float
Kernel width in respective dimension.
"""
markers = args_markers.markers
Np = args_markers.Np
out = 0.0
for neigh in range(27):
box_to_search = neighbours[loc_box, neigh]
c = 0
# loop over all particles in a box
while boxes[box_to_search, c] != -1:
p = boxes[box_to_search, c]
c += 1
if not holes[p]:
r1 = distance(eta1, markers[0, p], periodic1)
r2 = distance(eta2, markers[1, p], periodic2)
r3 = distance(eta3, markers[2, p], periodic3)
out += markers[index, p] * sph_smoothing_kernels.smoothing_kernel(kernel_type, r1, r2, r3, h1, h2, h3)
return out
####################
# naive evaluation #
####################
########################
# box-based evaluation #
########################
