锈蚀过程中的光滑粒子流体力学

Question

所以我同时学习了锈蚀和平滑的粒子流体力学。

我一直在使用来自AMD的这段视频作为参考和提供的平滑内核。

在这个阶段，我没有用体素来优化模拟，也没有增加硬件加速，所以它非常慢。我也没有添加一个可视化步骤来检查它是否真的有效。

此时，它编译并运行。因此，我感兴趣的是，从纯粹的生锈代码的角度来看，这段代码看上去是否可以理解。如果你想回顾一下实际的SPH实现，我不会说不.

//use std::collections::HashMap;
use std::f32::consts::PI;
use cgmath::{InnerSpace, Vector3, Vector4, Zero};
use rand::Rng;

const NUM_PARTICLES: u32 = 1000; // 1.2 Litres
const BUFFER_LEN: usize = NUM_PARTICLES as usize;
const BOX_DIMENSIONS_M: Vector3<f32> = Vector3::new(0.30, 0.10, 0.10);
const SMOOTHING_RADIUS: f32 = NUM_PARTICLES as f32 / (0.3 * 0.1 * 0.1 * 60.0);
const TIME_STEP_SECONDS: f32 = 0.01;
const GRAVITY: Vector4<f32> = Vector4::new(0.0, -9.8, 0.0, 0.0);
const SIM_LENGTH_SECONDS: f32 = 10.0;
const PARTICLE_MASS_KG: f32 = 0.001; //1mL of water
const FLUID_CONST: f32 = 2.2; // Possibly needs a x10^9
const VISCOUS_CONST: f32 = 0.0001;
const RHO_ZERO: f32 = 1000.0;

//TODO - Voxelise

fn main() {
    let mut sim_data: SimulationData = setup();
    let mut current_time: f32 = 0.0;
    while current_time < SIM_LENGTH_SECONDS {
        sim_step(&mut sim_data);
        current_time += TIME_STEP_SECONDS;
        println!("Time: {}s", current_time);
    }
}

fn sim_step(sim_data: &mut SimulationData) {

    let densities = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, _)| calculate_density_at_point(&sim_data.simulation_space, i))
        .collect::<Vec<f32>>();

    let pressures = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, _)| calculate_pressure_at_point(densities[i]))
        .collect::<Vec<f32>>();

    let pressure_grad_terms = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, _)| calculate_pressure_grad_term_at_point(&sim_data.simulation_space, &pressures, &densities, i))
        .collect::<Vec<Vector4<f32>>>();

    let viscosity_terms = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, _)| calculate_viscosity_term_at_point(&sim_data.simulation_space, &densities, i))
        .collect::<Vec<Vector4<f32>>>();

    let accelerations = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, _)| calculate_acceleration_at_point(&pressure_grad_terms, &viscosity_terms, i))
        .collect::<Vec<Vector4<f32>>>();

    let new_velocities = sim_data.simulation_space.velocities.iter().enumerate()
        .map(|(i, v)| v + accelerations[i] * TIME_STEP_SECONDS)
        .collect::<Vec<Vector4<f32>>>();

    let new_positions = sim_data.simulation_space.positions.iter().enumerate()
        .map(|(i, p)| p + new_velocities[i] * TIME_STEP_SECONDS)
        .collect::<Vec<Vector4<f32>>>();

    //TODO - Check for collisions with walls and floor

    // This will make rustaceans cry
    sim_data.simulation_space.velocities = new_velocities.try_into().unwrap();
    sim_data.simulation_space.positions = new_positions.try_into().unwrap();
}

fn calculate_acceleration_at_point(pressure_terms: &Vec<Vector4<f32>>, viscosity_terms: &Vec<Vector4<f32>>, i: usize) -> Vector4<f32> {
    let acceleration: Vector4<f32> = GRAVITY + pressure_terms[i] + viscosity_terms[i];
    return acceleration;
}

fn calculate_viscosity_term_at_point(sim_space: &SimulationSpace, densities: &Vec<f32>, i: usize) -> Vector4<f32> {
    let viscosity_term: Vector4<f32> = sim_space.positions.iter()
        .filter(|&x| is_in_interaction_radius_and_not_self(x, sim_space.positions[i]))
        .enumerate()// exclude self
        .fold(Vector4::zero(), |acc: Vector4<f32>, (j, _)| {
            acc + (VISCOUS_CONST / densities[j]) * PARTICLE_MASS_KG * (sim_space.velocities[j] - sim_space.velocities[i]) / densities[j] * laplacian_smooth(sim_space.positions[i], sim_space.positions[j])
        });
    return viscosity_term;
}

fn calculate_pressure_at_point(density: f32) -> f32 {
    let pressure_at_point: f32 = FLUID_CONST * (density - RHO_ZERO);
    return pressure_at_point;
}

fn calculate_pressure_grad_term_at_point(sim_space: &SimulationSpace, pressures: &Vec<f32>, densities: &Vec<f32>, i: usize) -> Vector4<f32> {
    let pressure_grad_term: Vector4<f32> = sim_space.positions.iter()
        .filter(|&x| is_in_interaction_radius_and_not_self(x, sim_space.positions[i]))
        .enumerate()// exclude self
        .fold(Vector4::zero(), |acc: Vector4<f32>, (j, _)| {
            acc + PARTICLE_MASS_KG * (pressures[i] / (densities[i].powi(2)) + pressures[j] / (densities[j].powi(2))) * grad_smooth(sim_space.positions[i], sim_space.positions[j])
        });
    return pressure_grad_term;
}

fn calculate_density_at_point(sim_space: &SimulationSpace, i: usize) -> f32 {
    let density = sim_space.positions.iter()
        .filter(|&x| is_in_interaction_radius_and_not_self(x, sim_space.positions[i]))
        .enumerate()// exclude self
        .fold(0.0, |acc: f32, (j, _)| acc + PARTICLE_MASS_KG * smooth(sim_space.positions[i], sim_space.positions[j]));

    return density;
}

fn is_in_interaction_radius_and_not_self(current: &Vector4<f32>, other: Vector4<f32>) -> bool {
    return *current != other && (current - other).magnitude() < SMOOTHING_RADIUS;
}

// There are some terms reused between smooth, grad_smooth and laplacian_smooth, this could be optimised
// Not to mention, some of these terms are constants...
fn smooth(current_position: Vector4<f32>, other_position: Vector4<f32>) -> f32
{
    return (315.0/(64.0*PI*(SMOOTHING_RADIUS.powi(9)))) * (SMOOTHING_RADIUS.powi(2) - (current_position - other_position).magnitude2()).powi(3);
}

fn grad_smooth(current_position: Vector4<f32>, other_position: Vector4<f32>) -> Vector4<f32>
{

    return (-45.0/(PI*(SMOOTHING_RADIUS.powi(6)))) * (SMOOTHING_RADIUS - (current_position - other_position).magnitude2()).powi(2) * ((current_position - other_position) / (current_position - other_position).magnitude2());
}

fn laplacian_smooth(current_position: Vector4<f32>, other_position: Vector4<f32>) -> f32
{
    return (45.0/(PI*(SMOOTHING_RADIUS.powi(6)))) * (SMOOTHING_RADIUS - (current_position - other_position).magnitude());
}

fn setup() -> SimulationData {
    let sim_space = SimulationSpace{
        positions: get_initial_positions(),
        velocities: get_initial_velocities(),
        accelerations: get_initial_accelerations()
    };

    let sim_data = SimulationData {
        simulation_space: sim_space,
        //voxel_pixel_map: HashMap::new()
    };

    return sim_data;
}

fn get_initial_positions() -> [Vector4<f32>; BUFFER_LEN] {
    let mut rng = rand::thread_rng();
    let mut positions: [Vector4<f32>; BUFFER_LEN] = [Vector4::zero(); BUFFER_LEN];
    for i in 0..BUFFER_LEN {
        positions[i] = Vector4::new(
            rng.gen_range(0.0..0.1),
            rng.gen_range(0.0..0.01),
            rng.gen_range(4.9..5.0),
            0.0
        );
    }
    return positions;
}

fn get_initial_velocities() -> [Vector4<f32>; BUFFER_LEN] {
    return [ Vector4::zero(); BUFFER_LEN];
}

fn get_initial_accelerations() -> [Vector4<f32>; BUFFER_LEN] {
    return [ Vector4::zero(); BUFFER_LEN];
}

struct SimulationData {
    simulation_space: SimulationSpace,
    //voxel_pixel_map: HashMap<u16, Vec<u32>>
}

struct SimulationSpace {
    positions: [Vector4<f32>; BUFFER_LEN as usize],
    velocities: [Vector4<f32>; BUFFER_LEN as usize],
    accelerations: [Vector4<f32>; BUFFER_LEN as usize]
}

Answer 1

一个重要的项目是:在--release模式下编译和运行。您提到它是缓慢的，但它只需要几秒钟的运行，我在发布模式，主要是因为所有的打印它做。

通常，在Rust中，我们会将setup函数作为SimulationData上的构造函数

impl SimulationData {
   fn new() -> SimulationData {
      ...
   }
}

我们通常也会让sim_step成为SimulationData上的一个方法

impl SimulationData {
    fn step(&mut self) {
        ...
    }
}

这句话似乎难以理解，而且可能是错误的。

const NUM_PARTICLES: u32 = 1000; // 1.2 Litres

这个常量可以用一个注释来解释关于神秘数字被乘以的原因：

const SMOOTHING_RADIUS: f32 = NUM_PARTICLES as f32 / (0.3 * 0.1 * 0.1 * 60.0);

你有很多块有这样的模式

let pressure_grad_term: Vector4<f32> = sim_space.positions.iter()
    .filter(|&x| is_in_interaction_radius_and_not_self(x, sim_space.positions[i]))
    .enumerate()// exclude self
    .fold(Vector4::zero(), |acc: Vector4<f32>, (j, _)| {
        acc + PARTICLE_MASS_KG * (pressures[i] / (densities[i].powi(2)) + pressures[j] / (densities[j].powi(2))) * grad_smooth(sim_space.positions[i], sim_space.positions[j])
    });

首先，这是错误的，因为你filter在你enumerate之前。这意味着所有索引都关闭，因为索引是在筛选器之后生成的。如果希望索引与数组匹配，则需要在筛选出项之前进行enumerate()。

其次，铁锈有一个sum方法，您可以使用，而不是一般折叠。

 let pressure_grad_term: Vector4<f32> = sim_space.positions.iter()
    .filter(|&x| is_in_interaction_radius_and_not_self(x, sim_space.positions[i]))
    .enumerate()// exclude self
    .map(|(j, _)| {
        PARTICLE_MASS_KG * (pressures[i] / (densities[i].powi(2)) + pressures[j] / (densities[j].powi(2))) * grad_smooth(sim_space.positions[i], sim_space.positions[j])
    })
    .sum()

但是，对于您在这里所做的编码类型，您可能想看看如何使用ndarray机箱。它提供了一个Array类型，它具有对数组进行操作的有用函数。在您的例子中，这似乎比转换为迭代器和返回更合适。例如，下面是calculate_density_at_point的重写版本

fn calculate_density_at_point(sim_space: &SimulationSpace, i: usize) -> f32 {
    sim_space.positions.fold(0.0, |acc, x| {
        if is_in_interaction_radius_and_not_self(x, sim_space.positions[i]) {
            acc + PARTICLE_MASS_KG * smooth(sim_space.positions[i], *x)
        } else {
            acc
        }
    })
}

您会注意到，它比迭代器版本简单得多。我们可以更进一步，定义一个函数来构建整个密度数组。

fn calculate_density_at_point(sim_space: &SimulationSpace) -> Array1<f32> {
    sim_space.positions.map(|&position| {
        sim_space.positions.fold(0.0, |acc, x| {
            if is_in_interaction_radius_and_not_self(x, position) {
                acc + PARTICLE_MASS_KG * smooth(position, *x)
            } else {
                acc
            }
        })
    })
}

更进一步，如果打开并行特性，则可以轻松地将计算并行化：

fn calculate_density_at_point(sim_space: &SimulationSpace) -> Array1<f32> {
    ndarray::Zip::from(&sim_space.positions).par_map_collect(|&position| {
        sim_space.positions.fold(0.0, |acc, x| {
            if is_in_interaction_radius_and_not_self(x, position) {
                acc + PARTICLE_MASS_KG * smooth(position, *x)
            } else {
                acc
            }
        })
    })
}