///|
pub(all) enum ShapeKind {
NGP
CIC
TSC
} derive(Eq, Debug, ToJson)
///|
pub(all) struct ShapeContribution {
index : Int
weight : Double
} derive(Debug, ToJson)
///|
pub(all) struct ShapeWeights {
contributions : Array[ShapeContribution]
total : Double
count : Int
} derive(Debug, ToJson)
///|
fn shape_kernel(distance : Double, kind : ShapeKind) -> Double {
let d = distance.abs()
match kind {
NGP => if d < 0.5 { 1.0 } else { 0.0 }
CIC => if d < 1.0 { 1.0 - d } else { 0.0 }
TSC =>
if d < 0.5 {
0.75 - d * d
} else if d < 1.5 {
0.5 * (1.5 - d) * (1.5 - d)
} else {
0.0
}
}
}
///|
pub fn shape_support(kind : ShapeKind) -> Int {
match kind {
NGP => 1
CIC => 2
TSC => 3
}
}
///|
pub fn shape_weights(
grid : Grid1D,
x : Double,
kind : ShapeKind,
) -> ShapeWeights {
let location = locate_periodic(grid, x)
let center = location.value
let radius = match kind {
NGP => 0
CIC => 1
TSC => 2
}
let candidates : Array[ShapeContribution] = []
for offset in -radius..<(radius + 1) {
let candidate = @math.floor(center).to_int() + offset
let weight = shape_kernel(center - candidate.to_double(), kind)
if weight > 0.0 {
candidates.push({ index: modulo_index(candidate, grid.cells), weight })
}
}
let raw_total = candidates.fold(init=0.0, fn(acc, item) { acc + item.weight })
let normalized = candidates.map(fn(item) {
{ index: item.index, weight: safe_ratio(item.weight, raw_total, 0.0) }
})
{
contributions: normalized,
total: normalized.fold(init=0.0, fn(acc, item) { acc + item.weight }),
count: normalized.length(),
}
}
///|
pub fn shape_weight(
grid : Grid1D,
x : Double,
index : Int,
kind : ShapeKind,
) -> Double {
let weights = shape_weights(grid, x, kind)
let mut result = 0.0
for item in weights.contributions {
if item.index == modulo_index(index, grid.cells) {
result = result + item.weight
}
}
result
}
///|
pub fn deposit_with_shape(
grid : Grid1D,
positions : ArrayView[Double],
values : ArrayView[Double],
kind : ShapeKind,
) -> Array[Double] {
let output = zeros(grid.cells)
let count = if positions.length() < values.length() {
positions.length()
} else {
values.length()
}
for i in 0.. Array[Double] {
let output = zeros(positions.length())
for i in 0.. Array[Double] {
let values = particles.map(fn(p) { p.charge * p.weight / grid.dx })
let positions = particles.map(fn(p) { p.x })
deposit_with_shape(grid, positions, values, kind)
}
///|
pub fn deposit_number_with_shape(
grid : Grid1D,
particles : ArrayView[Particle],
kind : ShapeKind,
) -> Array[Double] {
let values = particles.map(fn(p) { p.weight / grid.dx })
let positions = particles.map(fn(p) { p.x })
deposit_with_shape(grid, positions, values, kind)
}
///|
pub fn shape_charge_error(
grid : Grid1D,
particles : ArrayView[Particle],
kind : ShapeKind,
) -> Double {
let rho = deposit_charge_with_shape(grid, particles, kind)
(total_charge(grid, rho) - total_particle_charge(particles)).abs()
}
///|
pub fn shape_moment(
grid : Grid1D,
positions : ArrayView[Double],
values : ArrayView[Double],
kind : ShapeKind,
order : Int,
) -> Double {
let deposited = deposit_with_shape(grid, positions, values, kind)
let mut result = 0.0
for i in 0.. Array[ShapeContribution] {
let total = weights.fold(init=0.0, fn(acc, item) { acc + item.weight })
weights.map(fn(item) {
{ index: item.index, weight: safe_ratio(item.weight, total, 0.0) }
})
}
///|
pub fn shape_name(kind : ShapeKind) -> String {
match kind {
NGP => "ngp"
CIC => "cic"
TSC => "tsc"
}
}
///|
pub fn shape_kind_from_index(index : Int) -> ShapeKind {
match modulo_index(index, 3) {
0 => NGP
1 => CIC
_ => TSC
}
}
///|
pub fn shape_blend(
grid : Grid1D,
x : Double,
first : ShapeKind,
second : ShapeKind,
fraction : Double,
) -> Array[Double] {
let a = shape_weights(grid, x, first)
let b = shape_weights(grid, x, second)
let output = zeros(grid.cells)
for item in a.contributions {
output[item.index] = output[item.index] + item.weight * (1.0 - fraction)
}
for item in b.contributions {
output[item.index] = output[item.index] + item.weight * fraction
}
output
}
///|
pub fn shape_reconstruct(
grid : Grid1D,
deposited : ArrayView[Double],
positions : ArrayView[Double],
kind : ShapeKind,
) -> Array[Double] {
gather_with_shape(grid, deposited, positions, kind)
}
///|
pub fn shape_partition_error(
grid : Grid1D,
x : Double,
kind : ShapeKind,
) -> Double {
let weights = shape_weights(grid, x, kind)
(weights.total - 1.0).abs()
}