|
Simple Molecular Dynamics of $N$ particles interacting via a linear spring potential.
This example creates $N$ particles within a two-dimensional square domain, with periodic boundary conditions.
There is a linear spring force $\mathbf{f}_{ij}$ between particles $i$ and $j$ with a rest separation of $r$ (constant for all particles), and a cutoff at $r$. That is, if $\mathbf{r}_i$ is the position of particle $i$ and $\mathbf{dx}_{ij}=\mathbf{r}_j-\mathbf{r}_j$, then
$$ \mathbf{f}_{ij} = \begin{cases} \frac{r-|\mathbf{dx}_{ij}|}{|\mathbf{dx}_{ij}|}\mathbf{dx}_{ij}, & \text{for } |\mathbf{dx}_{ij}|<r \\ 0 & \text{otherwise}. \end{cases} $$
We wish to use a leap frog integrator to evolve positions $\mathbf{r}_i$ using velocities $\mathbf{v}_i$ and accelerations $\mathbf{a}_i = \sum_j \mathbf{f}_{ij}$. This gives the following update equations for each timestep $n$
\begin{align*} \mathbf{v}^{n+1}_i &= \mathbf{v}^n_i + \frac{dt}{m_i} \sum_j \mathbf{f}^n_{ij} \\ \mathbf{r}^{n+1}_i &= \mathbf{r}^n_i + dt, \mathbf{v}^{n+1}_i. \end{align*}
We implement this in Aboria using the code given below. Firstly we create the particle set data structure and add particles, ensuring that we have an initial condition where all the spring forces are $\mathbf{f}_{ij}=0$. Then we start the timestep loop, using our update equations given above.
#include <boost/math/constants/constants.hpp> #include <math.h> #include <random> #include "Aboria.h" using namespace Aboria; int main() { const double PI = boost::math::constants::pi<double>(); // Create a 2d particle container with one additional variable // "velocity", represented by a 2d double vector ABORIA_VARIABLE(velocity, vdouble2, "velocity") typedef Particles<std::tuple<velocity>,2> container_type; typedef typename container_type::position position; // set parameters for the MD simulation const int timesteps = 3000; const int nout = 200; const int timesteps_per_out = timesteps / nout; const double L = 31.0 / 1000.0; const int N = 30; const double diameter = 0.0022; const double k = 1.0e01; const double dens = 1160.0; const double mass = PI * std::pow(0.5 * diameter, 2) * dens; const double reduced_mass = 0.5 * mass; const double dt = (1.0 / 50.0) * PI / std::sqrt(k / reduced_mass); // create N particles container_type particles(N); // create symbols and labels in order to use the expression API Symbol<position> p; Symbol<velocity> v; Symbol<id> id_; Uniform uniform; VectorSymbolic<double, 2> vector; Label<0, container_type> a(particles); Label<1, container_type> b(particles); // dx is a symbol representing the difference in positions of // particle a and b. auto dx = create_dx(a, b); // sum is a symbolic function that sums a sequence of 2d vectors AccumulateWithinDistance<std::plus<vdouble2>> sum(diameter); // initialise particles to a uniform distribution p[a] = L * vector(uniform[a], uniform[a]); // zero initial velocity v[a] = vector(0, 0); // initiate neighbour search on a periodic 2d domain of side length // L set average number of particles per cell to 1 particles.init_neighbour_search(vdouble2(0, 0), vdouble2(L, L), vbool2(true, true)); // perform MD timestepping for (int io = 0; io < nout; ++io) { // on every i/o step write particle container to a vtk // unstructured grid file std::cout << "." << std::flush; #ifdef HAVE_VTK vtkWriteGrid("particles", io, particles.get_grid(true)); #endif for (int i = 0; i < timesteps_per_out; i++) { // leap frog integrator v[a] += dt * // spring force between particles sum(b, if_else(id_[a] != id_[b], -k * (diameter / norm(dx) - 1.0), 0.0) * dx) / mass; // update positions p[a] += dt * v[a]; } } std::cout << std::endl; }
This example creates two particle sets, the first representing a set of point particles undergoing brownian motion, the second representing a set of spheres with a given radius, different for each sphere.
The point particles diffuse around the three-dimensional, periodic domain, and reflect off the spheres whenever they encounter them.
Let indicies $i$ and $j$ refer to a pair of point particles, and $a$ and $b$ refer to a pair of spheres. Let $\mathbf{p}_i$ be the position of particle $i$, and let $\mathbf{dx}_{ij}$ refer to the shortest difference between the positions of particles $i$ and $j$. Let $r_b$ be the radius of sphere $b$. Then the update equations used to evolve the system are given by
$$ \mathbf{p}_i = \mathbf{p}_i + \sqrt{2\, D\, dt}\, \mathbf{N} + \sum_b \begin{cases} -2\left(\frac{r_b}{||\mathbf{dx}_{ib}||} - 1\right) \mathbf{dx}_{ib}, & \text{for } ||\mathbf{dx}_{ib}||<r_b \0 & \text{otherwise}. \end{cases} $$
where $D$ is the diffusion constant, $dt$ is the timestep and $\mathbf{N}$ is a three-dimensional vector containing random samples from a normal distribution
#include "Aboria.h" #include <random> using namespace Aboria; int main() { ABORIA_VARIABLE(radius, double, "radius") typedef Particles<std::tuple<radius>,3,std::vector> spheres_type; typedef Particles<std::tuple<>,3,std::vector> points_type; typedef position_d<3> position; spheres_type spheres; const double L = 10.0; const double D = 1.0; const double dt = 0.01; const double timesteps = 500; spheres.push_back(vdouble3(0, 0, 0)); get<radius>(spheres[0]) = 1.0; spheres.push_back(vdouble3(5, 0, 0)); get<radius>(spheres[1]) = 2.0; spheres.push_back(vdouble3(0, -5, 0)); get<radius>(spheres[2]) = 1.5; spheres.push_back(vdouble3(0, 0, 5)); get<radius>(spheres[3]) = 1.0; points_type points; std::uniform_real_distribution<double> uni(-L + L / 5, L - L / 5); std::default_random_engine generator; for (int i = 0; i < 1000; ++i) { points.push_back( vdouble3(uni(generator), uni(generator), uni(generator))); } points.init_neighbour_search(vdouble3(-L, -L, -L), vdouble3(L, L, L), vbool3(true, true, true)); spheres.init_neighbour_search(vdouble3(-L, -L, -L), vdouble3(L, L, L), vbool3(false, false, false)); Symbol<position> p; Symbol<radius> r; Symbol<alive> alive_; Label<0, spheres_type> a(spheres); Label<1, spheres_type> b(spheres); Label<0, points_type> i(points); Label<1, points_type> j(points); auto dx = create_dx(i, b); Normal N; VectorSymbolic<double, 3> vector; AccumulateWithinDistance<std::bit_or<bool>> any(4); AccumulateWithinDistance<std::plus<vdouble3>> sum(4); int count_before = 0; for (auto point : points) { if ((get<position>(point) - get<position>(spheres[0])).norm() < get<radius>(spheres[0])) { count_before++; } } /* * Kill any points within spheres */ alive_[i] = !any(b, if_else(norm(dx) < r[b], true, false)); int count_after = 0; for (auto point : points) { if ((get<position>(point) - get<position>(spheres[0])).norm() < get<radius>(spheres[0])) { count_after++; } } std::cout << " found " << count_before << " before and " << count_after << " after" << std::endl; /* * Diffusion step for points and reflect off spheres */ for (int ts = 1; ts < timesteps; ++ts) { if (ts % 10 == 0) { #ifdef HAVE_VTK vtkWriteGrid("bd", ts / 10, points.get_grid(true)); #endif std::cout << "." << std::flush; } p[i] += std::sqrt(2 * D * dt) * vector(N[i], N[i], N[i]); p[i] += sum(b, if_else(norm(dx) < r[b] && norm(dx) != 0, -2 * (r[b] / norm(dx) - 1), 0) * dx); } std::cout << std::endl; }
An Discrete Element Model example, simulating a polydisperse set of spherical particles falling onto an surface.
We wish to describe a set of particles in 3D cuboid, which is periodic in the Cartesian $x$, and $y$ directions. There is a linear spring force $\mathbf{f}_{ij}$ plus dissipation term between particles $i$ and $j$ with a rest separation at $s_i + s_j$ and a cutoff at $s_i + s_j$. That is, if $\mathbf{r}_i$ is the position of particle $i$ and $\mathbf{dx}_{ij}=\mathbf{r}_j-\mathbf{r}_j$, then
$$ \mathbf{f}_{ij} = \begin{cases} \frac{s_i+s_j-|\mathbf{dx}_{ij}|}{|\mathbf{dx}_{ij}|}\mathbf{dx}_{ij} + \gamma(\mathbf{v}_j-\mathbf{v}_i), & \text{for } |\mathbf{dx}_{ij}|<s_i+s_j \\ 0 & \text{otherwise}. \end{cases} $$
We wish to use a leap frog integrator to evolve positions $\mathbf{r}_i$ using velocities $\mathbf{v}_i$ and accelerations $\mathbf{a}_i = \frac{1}{m_i} \sum_j \mathbf{f}_{ij} - \mathbf{g}$. This gives the following update equations for each timestep $n$
\begin{align*} \mathbf{v}^{n+1}_i &= \mathbf{v}^n_i + \frac{dt}{m_i} \sum_j \mathbf{f}^n_{ij} \\ \mathbf{r}^{n+1}_i &= \mathbf{r}^n_i + dt, \mathbf{v}^{n+1}_i. \end{align*}
This figure above shows the simulation domain. As well as periodic in $x$, and $y$, we wish to add a soft surface at $z=0$ that also interacts with the particles using the same linear spring force. The domain is left open at $z=h$. We wish to initialise $N$ particles within this domain, and ensure that they are non-interacting at $t=0$.
#include <random> #include "Aboria.h" using namespace Aboria; #include <boost/math/constants/constants.hpp> int main() { const double PI = boost::math::constants::pi<double>(); /* * create variable types */ ABORIA_VARIABLE(radius, double, "radius") ABORIA_VARIABLE(mass, double, "mass") ABORIA_VARIABLE(velocity, vdouble3, "velocity") ABORIA_VARIABLE(acceleration, vdouble3, "acceleration") /* * create particle container */ typedef Particles<std::tuple<radius, velocity, acceleration, mass>> dem_type; typedef position_d<3> position; dem_type dem; /* * simulation parameters */ const int timesteps = 3000; const int nout = 200; const int timesteps_per_out = timesteps / nout; const double L = 31.0 / 1000.0; const int N = 30; /* * dem parameters */ const double dem_diameter = 0.0022; const double dem_gamma = 0.0004; const double dem_k = 1.0e01; const double dem_dens = 1160.0; const double dem_mass_min = (1.0 / 6.0) * PI * std::pow(0.5 * dem_diameter, 3) * dem_dens; const double dem_mass_max = (1.0 / 6.0) * PI * std::pow(dem_diameter, 3) * dem_dens; const double dem_min_reduced_mass = dem_mass_min * dem_mass_max / (dem_mass_min + dem_mass_max); const double dt = (1.0 / 50.0) * PI / sqrt(dem_k / dem_min_reduced_mass - std::pow(0.5 * dem_gamma / dem_min_reduced_mass, 2)); /* * initialise neighbour search with 3d cuboid domain, * periodic in x and y, set cell size so that there is * an average of 2 particles within each cell */ dem.init_neighbour_search(vdouble3(0, 0, -dem_diameter), vdouble3(L / 3, L / 3, L + dem_diameter), vbool3(true, true, false), 2); /* * create N random, non-overlaping particles with * random radius between 0.5*dem_diameter and dem_diameter */ std::uniform_real_distribution<double> uni(0, 1); std::default_random_engine generator; for (int i = 0; i < N; ++i) { bool free_position = false; dem_type::value_type p; get<radius>(p) = 0.25 * (uni(generator) + 1) * dem_diameter; while (free_position == false) { get<position>(p) = vdouble3(uni(generator) * L / 3, uni(generator) * L / 3, uni(generator) * (L - dem_diameter) + dem_diameter / 2); free_position = true; /* * loop over all neighbouring particles within "dem_diameter" distance */ for (auto j = euclidean_search(dem.get_query(), get<position>(p), dem_diameter); j != false; ++j) { if (j.dx().norm() < get<radius>(*j) + get<radius>(p)) { free_position = false; break; } } } dem.push_back(p); } /* * create symbols and labels in order to use the expression API */ Symbol<position> p; Symbol<radius> r; Symbol<mass> m; Symbol<velocity> v; Symbol<acceleration> dvdt; Symbol<id> id_; Label<0, dem_type> a(dem); Label<1, dem_type> b(dem); /* * dx is a symbol representing the difference in positions of * particle a and b. */ auto dx = create_dx(a, b); /* * sum is a symbolic function that sums a sequence of 3d vectors * over neighbouring particles within "dem_diameter" */ AccumulateWithinDistance<std::plus<vdouble3>> sum(dem_diameter); /* * vector creates a new double vector of dimension 3 */ VectorSymbolic<double, 3> vector; /* * perform MD timestepping */ v[a] = vdouble3::Zero(); m[a] = (1.0 / 6.0) * PI * 8 * r[a] * r[a] * r[a] * dem_dens; for (int io = 0; io < nout; ++io) { std::cout << "." << std::flush; #ifdef HAVE_VTK vtkWriteGrid("dem", io, dem.get_grid(true), {{"time", io / 2.0}}); #endif for (int i = 0; i < timesteps_per_out; i++) { p[a] += v[a] * dt; dvdt[a] = ( // spring force between dem particles sum(b, if_else(norm(dx) < r[a] + r[b] && id_[a] != id_[b], -dem_k * ((r[a] + r[b]) / norm(dx) - 1) * dx + dem_gamma * (v[b] - v[a]), vdouble3::Zero())) // spring force between particles and bottom wall + if_else(r[a] - p[a][2] > 0, dem_k * (r[a] - p[a][2]), 0.0) * vdouble3(0, 0, 1) // gravity force + vdouble3(0, 0, -9.81) * m[a] ) / m[a]; v[a] += dvdt[a] * dt; } } std::cout << std::endl; }
An Smoothed Particle Hydrodynamics example, simulating a 2D water column over a no-slip boundary. The x and y directions are periodic.
Simulation of a water column in hydrostatic equilibrium. SPH discretises the Navier-Stokes equations using radial interpolation kernels defined over a given particle set. See the following papers for more details than can be described here:
SPH is based on the idea of kernel interpolation. A fluid variable $A(\mathbf{r})$ (such as velocity or density) is interpolated using a kernel $W$, which depends on the smoothing length variable $h$.
\begin{equation}\label{Eq:integralInterpolant} A(\mathbf{r}) = \int A(\mathbf{r'})W(\mathbf{r}-\mathbf{r'},h)d\mathbf{r'}. \end{equation}
where $m_b$ and $\rho_b$ are the mass and density of particle $b$.
Here we have used the fifth-order Wendland kernel for $W$, which has the form
\begin{eqnarray} W(q) = \frac{\beta}{h^d} \begin{cases} (2-q)^4(1+2q) & \text{for } 0 \leq q \le 2, \ 0 & \text{for } q > 2. \end{cases} \end{eqnarray}
where $q = ||\mathbf{r}-\mathbf{r'}||/h$.
The rate of change of density, or continuity equation, is given by
\begin{equation} \label{Eq:changeInDensity} \frac{D\rho_a}{Dt} = \frac{1}{\Omega_a}\sum_b
m_b \mathbf{v}_{ab} \cdot
abla_a W_{ab}, \end{equation}
where $\mathbf{v}_{ab}=\mathbf{v}_a-\mathbf{v}_b$. The correction term $\Omega_a$ due to variable $h$ is given by
\begin{equation} \Omega_a = 1 - \frac{\partial h_a}{\partial \rho_a} \sum_b m_b \frac{\partial W_{ab}(h_a)}{\partial h_a}. \end{equation}
The equation of state used is the standard quasi-compressible form
\begin{equation} P = B \left ( \left ( \frac{\rho}{\rho_0} \right )^\gamma - 1 \right ), \end{equation}
where $\gamma = 7$ is a typical value and $\rho_0$ is a reference density that is normally set to the density of the fluid.
The SPH momentum equation is given as below. Viscosity is included by adding a viscous term $\Pi$
\begin{equation}\label{Eq:sphJustPressureForce} \frac{d\mathbf{v}_a}{dt}
= -\sum_b m_b \left [ \left ( \frac{P_a}{\Omega_a \rho_a^2} + \Pi_{ab} \right
)
abla_a W_{ab}(h_a) + \left ( \frac{P_b}{\Omega_b \rho_b^2} + \Pi_{ab}
\right )
abla_a W_{ab}(h_b) \right ], \end{equation}
The SPH literature contains many different forms for $\Pi$. We have used the term
\begin{equation}\label{Eq:monaghansViscousTerm} \Pi_{ab} = - \alpha \frac{v_{sig} (\mathbf{v}_{ab} \cdot \mathbf{r}_{ab} )}{2 \overline{\rho}_{ab} |\mathbf{r}_{ab}|}, \end{equation}
where $v_{sig} = 2(c_s + |\mathbf{v}_{ab} \cdot \mathbf{r}_{ab}| / |\mathbf{r}_{ab}| )$ is a signal velocity that represents the speed at which information propagates between the particles.
The particle's position and velocity were integrated using the Leapfrog second order method, which is also reversible in time in the absence of viscosity. To preserve the reversibility of the simulation, $d\rho/dt$ was calculated using the particle's position and velocity at the end of the timestep, rather than the middle as is commonly done. The full integration scheme is given by
\begin{align} \mathbf{r}^{\frac{1}{2}} &= \mathbf{r}^{0} + \frac{\delta t}{2} \mathbf{v}^{0}, \\ \mathbf{v}^{\frac{1}{2}} &= \mathbf{v}^{0} + \frac{\delta t}{2} F(\mathbf{r}^{-\frac{1}{2}},\mathbf{v}^{-\frac{1}{2}},\rho^{-\frac{1}{2}}), \\ \rho^{\frac{1}{2}} &= \rho^{0} + \frac{\delta t}{2} D(\mathbf{r}^0,\mathbf{v}^0), \label{Eq:timestepDensity1} \\ \mathbf{v}^{1} &= \mathbf{v}^{0} + \delta t F(\mathbf{r}^{\frac{1}{2}},\mathbf{v}^{\frac{1}{2}},\rho^{\frac{1}{2}}), \\ \mathbf{r}^{1} &= \mathbf{r}^{\frac{1}{2}} + \frac{\delta t}{2} \mathbf{v}^{1}, \\ \rho^{1} &= \rho^{\frac{1}{2}} + \frac{\delta t}{2} D(\mathbf{r}^1,\mathbf{v}^1), \label{Eq:timestepDensity2} \end{align}
where $\mathbf{r}^0$, $\mathbf{r}^{1/2}$ and $\mathbf{r}^1$ is $\mathbf{r}$ at the start, mid-point and end of the timestep respectively. The functions $F$ and $D$ are respectivly the momentum and continuity equations given above. The timestep $\delta t$ is bounded by the standard Courant condition
\begin{equation} \delta t_1 \le \min_a \left ( 0.8 \frac{h_a}{v_{sig}} \right ), \end{equation}
where the minimum is taken over all the particles.
#include "Aboria.h" using namespace Aboria; #include <boost/math/constants/constants.hpp> const double PI = boost::math::constants::pi<double>(); const double WCON_WENDLAND = 21.0 / (256.0 * PI); const unsigned int NDIM = 3; /* * Note that we are using standard C++ function objects to implement * the kernel W and its gradient F. We can wrap these as lazy Aboria * functions that can be used within the symbolic Level 3 API. */ struct F_fun { typedef double result_type; double operator()(const double r, const double h) const { if (r == 0) return 0; const double q = r / h; if (q <= 2.0) { return (1 / std::pow(h, NDIM + 2)) * WCON_WENDLAND * (-4 * std::pow(2 - q, 3) * (1 + 2 * q) + 2 * std::pow(2 - q, 4)) / q; } else { return 0.0; } } }; struct W_fun { typedef double result_type; double operator()(const double r, const double h) const { const double q = r / h; if (q <= 2.0) { return (1 / std::pow(h, NDIM)) * WCON_WENDLAND * std::pow(2.0 - q, 4) * (1.0 + 2.0 * q); } else { return 0.0; } } }; ABORIA_BINARY_FUNCTION(F, F_fun, SymbolicDomain); ABORIA_BINARY_FUNCTION(W, W_fun, SymbolicDomain); int main() { ABORIA_VARIABLE(kernel_radius, double, "kernel radius"); ABORIA_VARIABLE(velocity, vdouble3, "velocity"); ABORIA_VARIABLE(velocity_tmp, vdouble3, "temp velocity"); ABORIA_VARIABLE(varh_omega, double, "varh omega"); ABORIA_VARIABLE(density, double, "density"); ABORIA_VARIABLE(total_force, vdouble3, "total force"); ABORIA_VARIABLE(is_fixed, uint8_t, "fixed boundary"); ABORIA_VARIABLE(pressure_div_density2, double, "pressure div density2"); typedef Particles< std::tuple<kernel_radius, velocity, velocity_tmp, varh_omega, density, total_force, is_fixed, pressure_div_density2>, 3> sph_type; typedef position_d<3> position; sph_type sph; Symbol<position> p; Symbol<id> id_; Symbol<velocity> v; Symbol<velocity_tmp> v0; Symbol<density> rho; Symbol<total_force> dvdt; Symbol<is_fixed> fixed; Symbol<varh_omega> omega; Symbol<kernel_radius> h; Symbol<pressure_div_density2> pdr2; Label<0, sph_type> a(sph); Label<1, sph_type> b(sph); const int timesteps = 200; const int nout = 10; const int timesteps_per_out = timesteps / nout; const double L = 31.0 / 1000.0; const int nx = 5; /* * sph parameters */ const double hfac = 1.5; const double visc = 8.9e-07; const double refd = 1000.0; const double dens = 1000.0; const double gamma = 7; const double VMAX = 2.0 * sqrt(2 * 9.81 * L); const double CSFAC = 10.0; const double spsound = CSFAC * VMAX; const double prb = std::pow(refd / dens, gamma - 1.0) * std::pow(spsound, 2) * refd / gamma; const double psep = L / nx; double dt = std::min(0.25 * hfac * psep / spsound, 0.125 * std::pow(hfac * psep, 2) / visc); const double mass = dens * std::pow(psep, NDIM); std::cout << "h = " << hfac * psep << " vmax = " << VMAX << std::endl; const double time_damping = dt * 500; double t = 0; const vdouble3 low(0, 0, -3.0 * psep); const vdouble3 high(L, L, L); const vbool3 periodic(true, true, false); for (int i = 0; i < nx; i++) { for (int j = 0; j < nx; j++) { for (int k = 0; k < nx + 3; k++) { typename sph_type::value_type p; get<position>(p) = low + vdouble3((i + 0.5) * psep, (j + 0.5) * psep, (k + 0.5) * psep); get<kernel_radius>(p) = hfac * psep; get<velocity>(p) = vdouble3(0, 0, 0); get<velocity_tmp>(p) = vdouble3(0, 0, 0); get<varh_omega>(p) = 1.0; get<density>(p) = dens; get<total_force>(p) = vdouble3(0, 0, 0); get<is_fixed>(p) = get<position>(p)[2] < 0; sph.push_back(p); } } } std::cout << "starting...." << std::endl; sph.init_neighbour_search(low, high, periodic); auto dx = create_dx(a, b); AccumulateWithinDistance<std::plus<vdouble3>> sum_vect(2 * hfac * psep); AccumulateWithinDistance<std::plus<double>> sum(2 * hfac * psep); Accumulate<Aboria::max<double>> max; max.set_init(0); Accumulate<Aboria::min<double>> min; min.set_init(L); VectorSymbolic<double, 3> vector; for (int i = 0; i < nout; ++i) { #ifdef HAVE_VTK vtkWriteGrid("sph", i, sph.get_grid(true)); #endif for (int k = 0; k < timesteps_per_out; ++k) { /* * 0 -> 1/2 step */ v[a] += if_else(fixed[a] == false, dt / 2, 0) * dvdt[a]; if (t < time_damping) v[a] *= 0.98; p[a] += dt / 2 * v[a]; /* * Calculate omega */ omega[a] = 1.0 - (mass / (rho[a] * NDIM)) * sum(b, if_else(norm(dx) < 2 * h[a], pow(norm(dx), 2) * F(norm(dx), h[a]) + NDIM * W(norm(dx), h[a]), 0)); /* * 1/2 -> 1 step */ /* * calculate change in density and kernel radius */ rho[a] += dt * (mass / omega[a]) * sum(b, if_else(norm(dx) < 2 * h[a], dot(v[b] - v[a], dx * F(norm(dx), h[a])), 0)); h[a] = pow(mass / rho[a], 1.0 / NDIM); /* * reset neighbour search for new kernel radius */ const double maxh = eval(max(a, h[a])); const double minh = eval(min(a, h[a])); sum_vect.set_max_distance(2 * maxh); sum.set_max_distance(2 * maxh); /* * advance velocity, position and calculate pdr2 */ v0[a] = v[a]; v[a] += if_else(fixed[a] == false, dt / 2.0, 0.0) * dvdt[a]; pdr2[a] = prb * (pow(rho[a] / refd, gamma) - 1.0) / pow(rho[a], 2); p[a] += dt / 2 * v0[a]; /* * acceleration due to gravity */ dvdt[a] = vdouble3(0, 0, -9.81); /* * acceleration on SPH calculation */ dvdt[a] += mass * sum_vect(b, if_else(norm(dx) < 2 * h[a], // pressure force ((1.0 / omega[a]) * pdr2[a] * F(norm(dx), h[a]) + (1.0 / omega[b]) * pdr2[b] * F(norm(dx), h[b])) * dx // viscous force + (v[a] - v[b]) * visc * (rho[a] + rho[b]) / (rho[a] * rho[b]) * F(norm(dx), h[a]) , vector(0, 0, 0))); /* * 1/2 -> 1 step for velocity */ v[a] = if_else(fixed[a] == false, v0[a] + dt / 2 * dvdt[a], vector(0, 0, 0)); t += dt; /* * sph timestep condition */ dt = std::min(0.25 * minh / spsound, 0.125 * std::pow(minh, 2) / visc); } std::cout << "iteration " << i << std::endl; } }