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