Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

1D support #249

Merged
merged 9 commits into from
Oct 31, 2023
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
184 changes: 76 additions & 108 deletions src/general/smoothing_kernels.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,67 @@
@doc raw"""
SmoothingKernel{NDIMS}()
LasNikas marked this conversation as resolved.
Show resolved Hide resolved

We provide three kernel functions, generated as the Fourier transform
```math
M_n(x, h) = \frac{1}{2 \pi}\int_{-\infty}^{\infty}
\left[ \frac{sin(k h /2)}{k h /2}\right]^n cos(k x) \mathrm{d}k \: .
```
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
| B-Spline | | compact support |
| -------- | :---------------------------------------- | ----------------- |
| $M_4$ | [`SchoenbergCubicSplineKernel`](@ref) | $[0, 2h]$ |
| $M_5$ | [`SchoenbergQuarticSplineKernel`](@ref) | $[0, 2.5h]$ |
| $M_6$ | [`SchoenbergQuinticSplineKernel`](@ref) | $[0, 3h]$ |
LasNikas marked this conversation as resolved.
Show resolved Hide resolved

In general, the Kernels are defined as
```math
W(\Vert r_a - r_b \Vert, h) = \frac{1}{h^d} w(q) \: ,
```
where $h$ is the smoothing length, $d$ refers to the number of spatial dimensions and
$q=\Vert r_a - r_b \Vert/h$.
LasNikas marked this conversation as resolved.
Show resolved Hide resolved

For an overview of Schoenberg cubic ($M_4$), quartic ($M_5$) and quintic ($M_6$) spline
kernels including normalization factors, see (Price, 2012).
For an analytic formula for higher order kernels, see (Monaghan, 1985).

LasNikas marked this conversation as resolved.
Show resolved Hide resolved
!!! note "Usage"
The kernel can be called as
```
TrixiParticles.kernel(::SmoothingKernel{NDIMS}, r, h)
```
The length of the compact support can be obtained as
```
TrixiParticles.compact_support(::SmoothingKernel{NDIMS}, h)
```

Note that ``r`` has to be a scalar, so in the context of SPH, the kernel
should be used as
```math
W(\Vert r_a - r_b \Vert, h).
```

The gradient required in SPH,
```math
\frac{\partial}{\partial r_a} W(\Vert r_a - r_b \Vert, h)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
```
can be called as
```
TrixiParticles.kernel_grad(kernel, pos_diff, distance, h)

```
where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$.

## References:
- Daniel J. Price. "Smoothed particle hydrodynamics and magnetohydrodynamics".
In: Journal of Computational Physics 231.3 (2012), pages 759-794.
[doi: 10.1016/j.jcp.2010.12.011](https://doi.org/10.1016/j.jcp.2010.12.011)
- Joseph J. Monaghan. "Particle methods for hydrodynamics".
In: Computer Physics Reports 3.2 (1985), pages 71–124.
[doi: 10.1016/0167-7977(85)90010-3](https://doi.org/10.1016/0167-7977(85)90010-3)
- Isaac J. Schoenberg. "Contributions to the problem of approximation of equidistant data by analytic functions.
Part B. On the problem of osculatory interpolation. A second class of analytic approximation formulae."
In: Quarterly of Applied Mathematics 4.2 (1946), pages 112–141.
[doi: 10.1090/QAM/16705](https://doi.org/10.1090/QAM/16705)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
"""
abstract type SmoothingKernel{NDIMS} end
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
@inline Base.ndims(::SmoothingKernel{NDIMS}) where {NDIMS} = NDIMS

Expand All @@ -21,51 +85,18 @@ end

Cubic spline kernel by Schoenberg (Schoenberg, 1946), given by
```math
W(r, h) = \frac{1}{h^d} w(r/h)
```
with
```math
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
w(q) = \sigma \begin{cases}
\frac{1}{4} (2 - q)^3 - (1 - q)^3 & \text{if } 0 \leq q < 1, \\
\frac{1}{4} (2 - q)^3 & \text{if } 1 \leq q < 2, \\
0 & \text{if } q \geq 2, \\
\end{cases}
```
where ``d`` is the number of dimensions and ``\sigma = 17/(7\pi)``
in two dimensions or ``\sigma = 1/\pi`` in three dimensions is a normalization factor.
where ``\sigma`` is a normalisation constant given by
$\sigma =[\frac{2}{3}, \frac{10}{7 \pi}, \frac{1}{\pi}]$ in $[1, 2, 3]$ dimensions.

This kernel function has a compact support of ``[0, 2h]``.

For an overview of Schoenberg cubic, quartic and quintic spline kernels
including normalization factors, see (Price, 2012).
For an analytic formula for higher order kernels, see (Monaghan, 1985).

!!! note "Usage"
The kernel can be called as `TrixiParticles.kernel(::SchoenbergCubicSplineKernel, r, h)`.
The length of the compact support can be obtained as
`TrixiParticles.compact_support(::SchoenbergCubicSplineKernel, h)`.

Note that ``r`` has to be a scalar, so in the context of SPH, the kernel
should be used as ``W(\Vert r_a - r_b \Vert, h)``.
The gradient required in SPH,
```math
\frac{\partial}{\partial r_a} W(\Vert r_a - r_b \Vert, h)
```
can be called as
`TrixiParticles.kernel_grad(kernel, pos_diff, distance, h)`,
where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$.

## References:
- Daniel J. Price. "Smoothed particle hydrodynamics and magnetohydrodynamics".
In: Journal of Computational Physics 231.3 (2012), pages 759-794.
[doi: 10.1016/j.jcp.2010.12.011](https://doi.org/10.1016/j.jcp.2010.12.011)
- Joseph J. Monaghan. "Particle methods for hydrodynamics".
In: Computer Physics Reports 3.2 (1985), pages 71–124.
[doi: 10.1016/0167-7977(85)90010-3](https://doi.org/10.1016/0167-7977(85)90010-3)
- Isaac J. Schoenberg. "Contributions to the problem of approximation of equidistant data by analytic functions.
Part B. On the problem of osculatory interpolation. A second class of analytic approximation formulae."
In: Quarterly of Applied Mathematics 4.2 (1946), pages 112–141.
[doi: 10.1090/QAM/16705](https://doi.org/10.1090/QAM/16705)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
For general information see [`SmoothingKernel`](@ref).
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
"""
struct SchoenbergCubicSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end

Expand Down Expand Up @@ -104,6 +135,7 @@ end

@inline compact_support(::SchoenbergCubicSplineKernel, h) = 2 * h

@inline normalization_factor(::SchoenbergCubicSplineKernel{1}, h) = 2 / (3 * h)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
@inline normalization_factor(::SchoenbergCubicSplineKernel{2}, h) = 10 / (7 * pi * h^2)
@inline normalization_factor(::SchoenbergCubicSplineKernel{3}, h) = 1 / (pi * h^3)

Expand All @@ -112,10 +144,6 @@ end

Quartic spline kernel by Schoenberg (Schoenberg, 1946), given by
```math
W(r, h) = \frac{1}{h^d} w(r/h)
```
with
```math
w(q) = \sigma \begin{cases}
\left(5/2 - q \right)^4 - 5\left(3/2 - q \right)^4
+ 10\left(1/2 - q \right)^4 & \text{if } 0 \leq q < \frac{1}{2}, \\
Expand All @@ -125,41 +153,12 @@ w(q) = \sigma \begin{cases}
0 & \text{if } q \geq \frac{5}{2},
\end{cases}
```
where ``d`` is the number of dimensions and ``\sigma = 96/(1199\pi)``
in two dimensions or ``\sigma = 1/(20\pi)`` in three dimensions is a normalization factor.
where ``\sigma`` is a normalisation constant given by
$\sigma =[\frac{1}{24}, \frac{96}{1199 \pi}, \frac{1}{20\pi}]$ in $[1, 2, 3]$ dimensions.

This kernel function has a compact support of ``[0, 2.5h]``.

For an overview of Schoenberg cubic, quartic and quintic spline kernels
including normalization factors, see (Price, 2012).
For an analytic formula for higher order kernels, see (Monaghan, 1985).

!!! note "Usage"
The kernel can be called as `TrixiParticles.kernel(::SchoenbergQuarticSplineKernel, r, h)`.
The length of the compact support can be obtained as
`TrixiParticles.compact_support(::SchoenbergQuarticSplineKernel, h)`.

Note that ``r`` has to be a scalar, so in the context of SPH, the kernel
should be used as ``W(\Vert r_a - r_b \Vert, h)``.
The gradient required in SPH,
```math
\frac{\partial}{\partial r_a} W(\Vert r_a - r_b \Vert, h)
```
can be called as
`TrixiParticles.kernel_grad(kernel, pos_diff, distance, h)`,
where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$.

## References:
- Daniel J. Price. "Smoothed particle hydrodynamics and magnetohydrodynamics".
In: Journal of Computational Physics 231.3 (2012), pages 759-794.
[doi: 10.1016/j.jcp.2010.12.011](https://doi.org/10.1016/j.jcp.2010.12.011)
- Joseph J. Monaghan. "Particle methods for hydrodynamics".
In: Computer Physics Reports 3.2 (1985), pages 71–124.
[doi: 10.1016/0167-7977(85)90010-3](https://doi.org/10.1016/0167-7977(85)90010-3)
- Isaac J. Schoenberg. "Contributions to the problem of approximation of equidistant data by analytic functions.
Part B. On the problem of osculatory interpolation. A second class of analytic approximation formulae."
In: Quarterly of Applied Mathematics 4.2 (1946), pages 112–141.
[doi: 10.1090/QAM/16705](https://doi.org/10.1090/QAM/16705)
For general information see [`SmoothingKernel`](@ref).
"""
struct SchoenbergQuarticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end

Expand Down Expand Up @@ -206,6 +205,7 @@ end

@inline compact_support(::SchoenbergQuarticSplineKernel, h) = 2.5 * h

@inline normalization_factor(::SchoenbergQuarticSplineKernel{1}, h) = 1 / (24 * h)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
@inline normalization_factor(::SchoenbergQuarticSplineKernel{2}, h) = 96 / (1199 * pi * h^2)
@inline normalization_factor(::SchoenbergQuarticSplineKernel{3}, h) = 1 / (20 * pi * h^3)

Expand All @@ -214,52 +214,19 @@ end

Quintic spline kernel by Schoenberg (Schoenberg, 1946), given by
```math
W(r, h) = \frac{1}{h^d} w(r/h)
```
with
```math
w(q) = \sigma \begin{cases}
(3 - q)^5 - 6(2 - q)^5 + 15(1 - q)^5 & \text{if } 0 \leq q < 1, \\
(3 - q)^5 - 6(2 - q)^5 & \text{if } 1 \leq q < 2, \\
(3 - q)^5 & \text{if } 2 \leq q < 3, \\
0 & \text{if } q \geq 3,
\end{cases}
```
where ``d`` is the number of dimensions and ``\sigma = 7/(478\pi)``
in two dimensions or ``\sigma = 1/(120\pi)`` in three dimensions is a normalization factor.
where ``\sigma`` is a normalisation constant given by
$\sigma =[\frac{1}{120}, \frac{7}{478 \pi}, \frac{1}{120\pi}]$ in $[1, 2, 3]$ dimensions.

This kernel function has a compact support of ``[0, 3h]``.

For an overview of Schoenberg cubic, quartic and quintic spline kernels
including normalization factors, see (Price, 2012).
For an analytic formula for higher order kernels, see (Monaghan, 1985).

!!! note "Usage"
The kernel can be called as `TrixiParticles.kernel(::SchoenbergQuinticSplineKernel, r, h)`.
The length of the compact support can be obtained as
`TrixiParticles.compact_support(::SchoenbergQuinticSplineKernel, h)`.

Note that ``r`` has to be a scalar, so in the context of SPH, the kernel
should be used as ``W(\Vert r_a - r_b \Vert, h)``.
The gradient required in SPH,
```math
\frac{\partial}{\partial r_a} W(\Vert r_a - r_b \Vert, h)
```
can be called as
`TrixiParticles.kernel_grad(kernel, pos_diff, distance, h)`,
where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$.

## References:
- Daniel J. Price. "Smoothed particle hydrodynamics and magnetohydrodynamics".
In: Journal of Computational Physics 231.3 (2012), pages 759-794.
[doi: 10.1016/j.jcp.2010.12.011](https://doi.org/10.1016/j.jcp.2010.12.011)
- Joseph J. Monaghan. "Particle methods for hydrodynamics".
In: Computer Physics Reports 3.2 (1985), pages 71–124.
[doi: 10.1016/0167-7977(85)90010-3](https://doi.org/10.1016/0167-7977(85)90010-3)
- Isaac J. Schoenberg. "Contributions to the problem of approximation of equidistant data by analytic functions.
Part B. On the problem of osculatory interpolation. A second class of analytic approximation formulae."
In: Quarterly of Applied Mathematics 4.2 (1946), pages 112–141.
[doi: 10.1090/QAM/16705](https://doi.org/10.1090/QAM/16705)
For general information see [`SmoothingKernel`](@ref).
"""
struct SchoenbergQuinticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end

Expand Down Expand Up @@ -306,5 +273,6 @@ end

@inline compact_support(::SchoenbergQuinticSplineKernel, h) = 3 * h

@inline normalization_factor(::SchoenbergQuinticSplineKernel{1}, h) = 1 / (120 * h)
LasNikas marked this conversation as resolved.
Show resolved Hide resolved
@inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (478 * pi * h^2)
@inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (120 * pi * h^3)
11 changes: 11 additions & 0 deletions src/neighborhood_search/grid_nhs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,17 @@ end
end
end

@inline function eachneighbor(coords, neighborhood_search::GridNeighborhoodSearch{1})
cell = cell_coords(coords, neighborhood_search)
x = cell[1]
# Generator of all neighboring cells to consider
neighboring_cells = ((x + i) for i in -1:1)

# Merge all lists of particles in the neighboring cells into one iterator
Iterators.flatten(particles_in_cell(cell, neighborhood_search)
for cell in neighboring_cells)
end

@inline function eachneighbor(coords, neighborhood_search::GridNeighborhoodSearch{2})
cell = cell_coords(coords, neighborhood_search)
x, y = cell
Expand Down
25 changes: 25 additions & 0 deletions src/setups/rectangular_shape.jl
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,24 @@ function rectangular_shape_coords(particle_spacing, n_particles_per_dimension,
return coordinates
end

# 1D
function initialize_rectangular!(coordinates, particle_spacing,
min_coordinates,
n_particles_per_dimension::NTuple{1}, loop_order)
n_particles_x = n_particles_per_dimension[1]
particle = 0

if loop_order === :x_first
for x in 1:n_particles_x
particle += 1
fill_coordinates!(coordinates, particle, min_coordinates, x, particle_spacing)
end

else
throw(ArgumentError("$loop_order is not a valid loop order. Possible values are :x_first"))
end
end

# 2D
function initialize_rectangular!(coordinates, particle_spacing,
min_coordinates,
Expand Down Expand Up @@ -141,6 +159,13 @@ function initialize_rectangular!(coordinates, particle_spacing,
end
end

@inline function fill_coordinates!(coordinates, particle,
min_coordinates, x, particle_spacing)
# The first particle starts at a distance `0.5particle_spacing` from `min_coordinates`
# in each dimension.
coordinates[1, particle] = min_coordinates[1] + (x - 0.5) * particle_spacing
end

@inline function fill_coordinates!(coordinates, particle,
min_coordinates, x, y, particle_spacing)
# The first particle starts at a distance `0.5particle_spacing` from `min_coordinates`
Expand Down