1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

#[cfg(feature = "parallel")]
use rayon::prelude::*;

use crate::geometry::ParticlesContacts;

use crate::math::{Real, Vector};
use crate::object::{Boundary, Fluid};
use crate::solver::NonPressureForce;
use crate::TimestepManager;

// See http://www.astro.lu.se/~david/teaching/SPH/notes/annurev.aa.30.090192.pdf
/// Implements artificial viscosity.
#[derive(Clone)]
pub struct ArtificialViscosity {
    /// The coefficient of the linear part of the viscosity.
    pub alpha: Real,
    /// The coefficient of the quadratic part of the viscosity.
    pub beta: Real,
    /// The speed of sound.
    pub speed_of_sound: Real,
    /// The fluid viscosity coefficient.
    pub fluid_viscosity_coefficient: Real,
    /// The viscosity coefficient when interacting with boundaries.
    pub boundary_viscosity_coefficient: Real,
}

impl ArtificialViscosity {
    /// Initializes the artificial viscosity with the given viscosity coefficients.
    pub fn new(fluid_viscosity_coefficient: Real, boundary_viscosity_coefficient: Real) -> Self {
        Self {
            alpha: na::one::<Real>(),
            beta: na::convert::<_, Real>(0.0),
            speed_of_sound: na::convert::<_, Real>(10.0),
            fluid_viscosity_coefficient,
            boundary_viscosity_coefficient,
        }
    }
}

impl NonPressureForce for ArtificialViscosity {
    fn solve(
        &mut self,
        _timestep: &TimestepManager,
        kernel_radius: Real,
        fluid_fluid_contacts: &ParticlesContacts,
        fluid_boundaries_contacts: &ParticlesContacts,
        fluid: &mut Fluid,
        boundaries: &[Boundary],
        densities: &[Real],
    ) {
        let fluid_viscosity_coefficient = self.fluid_viscosity_coefficient;
        let boundary_viscosity_coefficient = self.boundary_viscosity_coefficient;
        let speed_of_sound = self.speed_of_sound;
        let alpha = self.alpha;
        let beta = self.beta;
        let density0 = fluid.density0;
        let volumes = &fluid.volumes;
        let positions = &fluid.positions;
        let velocities = &fluid.velocities;
        let _0_5: Real = na::convert::<_, Real>(0.5);

        par_iter_mut!(fluid.accelerations)
            .enumerate()
            .for_each(|(i, acceleration)| {
                let mut fluid_acc = Vector::zeros();
                let mut boundary_acc = Vector::zeros();

                if self.fluid_viscosity_coefficient != na::zero::<Real>() {
                    for c in fluid_fluid_contacts
                        .particle_contacts(i)
                        .read()
                        .unwrap()
                        .iter()
                    {
                        if c.i_model == c.j_model {
                            let r_ij = positions[c.i] - positions[c.j];
                            let v_ij = velocities[c.i] - velocities[c.j];
                            let vr = r_ij.dot(&v_ij);

                            if vr < na::zero::<Real>() {
                                let density_average = (densities[c.i] + densities[c.j]) * _0_5;
                                let eta2 =
                                    kernel_radius * kernel_radius * na::convert::<_, Real>(0.01);
                                let mu_ij = kernel_radius * vr / (r_ij.norm_squared() + eta2);

                                fluid_acc += c.gradient
                                    * (fluid_viscosity_coefficient
                                        * (speed_of_sound * alpha * mu_ij - beta * mu_ij * mu_ij)
                                        * (volumes[c.j] * density0 / density_average));
                            }
                        }
                    }
                }

                if self.boundary_viscosity_coefficient != na::zero::<Real>() {
                    for c in fluid_boundaries_contacts
                        .particle_contacts(i)
                        .read()
                        .unwrap()
                        .iter()
                    {
                        let r_ij = positions[c.i] - boundaries[c.j_model].positions[c.j];
                        let v_ij = velocities[c.i] - boundaries[c.j_model].velocities[c.j];
                        let vr = r_ij.dot(&v_ij);

                        if vr < na::zero::<Real>() {
                            let density_average = densities[c.i];
                            let eta2 = kernel_radius * kernel_radius * na::convert::<_, Real>(0.01);
                            let mu_ij = kernel_radius * vr / (r_ij.norm_squared() + eta2);

                            boundary_acc += c.gradient
                                * (boundary_viscosity_coefficient
                                    * (speed_of_sound * alpha * mu_ij - beta * mu_ij * mu_ij)
                                    * (boundaries[c.j_model].volumes[c.j] * density0
                                        / density_average));
                            let mi = volumes[c.i] * density0;
                            boundaries[c.j_model].apply_force(c.j, boundary_acc * -mi);
                        }
                    }
                }

                *acceleration += fluid_acc + boundary_acc;
            })
    }

    fn apply_permutation(&mut self, _: &[usize]) {}
}