1 #include <assert.h>
2 #include <math.h>
3 
4 #include "kernels.h"
5 
6 #ifndef M_PI
7 #define M_PI 3.14159265358979323846
8 #endif
9 
compute_density(int n,float * restrict rho,float * restrict pos_x,float * restrict pos_y,float h,float mass)10 void compute_density(int n, float *restrict rho, float *restrict pos_x, float *restrict pos_y, float h, float mass)
11 {
12     float h_squared = h * h;
13     float h_pow_8 = h_squared * h_squared * h_squared * h_squared;
14     float C = 4 * mass / M_PI / h_pow_8;
15 
16     for (int i = 0; i < n; ++i) {
17         rho[i] = 4 * mass / M_PI / h_squared;
18     }
19 
20     for (int i = 0; i < n; ++i) {
21         for (int j = i + 1; j < n; ++j) {
22             float delta_x = pos_x[i] - pos_x[j];
23             float delta_y = pos_y[i] - pos_y[j];
24             float dist_squared = delta_x * delta_x + delta_y * delta_y;
25             float overlap = h_squared - dist_squared;
26 
27             if (overlap > 0) {
28                 float rho_ij = C * overlap * overlap * overlap;
29                 rho[i] += rho_ij;
30                 rho[j] += rho_ij;
31             }
32         }
33     }
34 }
35 
compute_accel(int n,float * restrict rho,float * restrict pos_x,float * restrict pos_y,float * restrict v_x,float * restrict v_y,float * restrict a_x,float * restrict a_y,float mass,float g,float h,float k,float rho0,float mu)36 void compute_accel(
37     int n,
38     float *restrict rho,
39     float *restrict pos_x,
40     float *restrict pos_y,
41     float *restrict v_x,
42     float *restrict v_y,
43     float *restrict a_x,
44     float *restrict a_y,
45     float mass,
46     float g,
47     float h,
48     float k,
49     float rho0,
50     float mu)
51 {
52     const float h_squared = h * h;
53     const float C_0 = mass / M_PI / (h_squared * h_squared);
54     const float C_p = 15 * k;
55     const float C_v = -40 * mu;
56 
57     // gravity:
58     for (int i = 0; i < n; ++i) {
59         a_x[i] = 0;
60         a_y[i] = -g;
61     }
62 
63     // Now compute interaction forces
64     for (int i = 0; i < n; ++i) {
65         for (int j = i + 1; j < n; ++j) {
66             float delta_x = pos_x[i] - pos_x[j];
67             float delta_y = pos_y[i] - pos_y[j];
68             float dist_squared = delta_x * delta_x + delta_y * delta_y;
69 
70             if (dist_squared < h_squared) {
71                 float q = sqrt(dist_squared) / h;
72                 float u = 1 - q;
73                 float w_0 = C_0 * u / rho[i] / rho[j];
74                 float w_p = w_0 * C_p * (rho[i] + rho[j] - 2 * rho0) * u / q;
75                 float w_v = w_0 * C_v;
76                 float delta_v_x = v_x[i] - v_x[j];
77                 float delta_v_y = v_y[i] - v_y[j];
78 
79                 a_x[i] += (w_p * delta_x + w_v * delta_v_x);
80                 a_y[i] += (w_p * delta_y + w_v * delta_v_y);
81                 a_x[j] -= (w_p * delta_x + w_v * delta_v_x);
82                 a_y[j] -= (w_p * delta_y + w_v * delta_v_y);
83             }
84         }
85     }
86 }
87 
damp_reflect(int which,float barrier,float * pos_x,float * pos_y,float * v_x,float * v_y)88 void damp_reflect(
89     int which,
90     float barrier,
91     float *pos_x,
92     float *pos_y,
93     float *v_x,
94     float *v_y)
95 {
96     float *v_which   = (which == 0) ? v_x   : v_y;
97     float *pos_which = (which == 0) ? pos_x : pos_y;
98 
99     // Coefficient of resitiution
100     const float DAMP = 0.75;
101     // Ignore degenerate cases
102     if (fabs(v_which[0]) <= 1e-3)
103         return;
104 
105     // Scale back the distance traveled based on time from collision
106     float tbounce = (pos_which[0] - barrier) / v_which[0];
107     pos_x[0] -= v_x[0]*(1-DAMP)*tbounce;
108     pos_y[0] -= v_y[0]*(1-DAMP)*tbounce;
109 
110     // Reflect the position and velocity
111     pos_which[0] = 2 * barrier - pos_which[0];
112     v_which[0]   = -v_which[0];
113 
114     // Damp the velocities
115     v_x[0] *= DAMP;
116     v_y[0] *= DAMP;
117 }
118 
reflect_bc(int n,float * restrict pos_x,float * restrict pos_y,float * restrict v_x,float * restrict v_y)119 void reflect_bc(
120     int n,
121     float *restrict pos_x,
122     float *restrict pos_y,
123     float *restrict v_x,
124     float *restrict v_y)
125 {
126     // Boundaries of the computational domain
127     const float XMIN = 0.0;
128     const float XMAX = 1.0;
129     const float YMIN = 0.0;
130     const float YMAX = 1.0;
131 
132     for (int i = 0; i < n; ++i) {
133         if (pos_x[i] < XMIN) {
134             damp_reflect(0, XMIN, pos_x + i, pos_y + i, v_x + i, v_y + i);
135         }
136         if (pos_x[i] > XMAX) {
137             damp_reflect(0, XMAX, pos_x + i, pos_y + i, v_x + i, v_y + i);
138         }
139         if (pos_y[i] < YMIN) {
140             damp_reflect(1, YMIN, pos_x + i, pos_y + i, v_x + i, v_y + i);
141         }
142         if (pos_y[i] > YMAX) {
143             damp_reflect(1, YMAX, pos_x + i, pos_y + i, v_x + i, v_y + i);
144         }
145     }
146 }
147 
leapfrog(int n,float * restrict pos_x,float * restrict pos_y,float * restrict v_x,float * restrict v_y,float * restrict a_x,float * restrict a_y,double dt)148 void leapfrog(
149     int n,
150     float *restrict pos_x,
151     float *restrict pos_y,
152     float *restrict v_x,
153     float *restrict v_y,
154     float *restrict a_x,
155     float *restrict a_y,
156     double dt)
157 {
158     for (int i = 0; i < n; ++i) {
159         v_x[i] += a_x[i] * dt;
160         v_y[i] += a_y[i] * dt;
161 
162         pos_x[i] += v_x[i] * dt;
163         pos_y[i] += v_y[i] * dt;
164     }
165 }
166