tsimulation.c - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tsimulation.c (5058B)
---
1 #ifndef __OpenBSD__
2 #define _POSIX_C_SOURCE 200809L
3 #endif
4 #include <unistd.h>
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <err.h>
8 #include <string.h>
9 #include <math.h>
10 #include "grain.h"
11 #include "grains.h"
12 #include "simulation.h"
13 #include "arrays.h"
14 #include "util.h"
15
16 struct simulation
17 sim_new(void)
18 {
19 struct simulation sim;
20 sim_defaults(&sim);
21 return sim;
22 }
23
24 void
25 sim_defaults(struct simulation *sim)
26 {
27 int d, ret;
28
29 ret = snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
30 if (ret < 0 || (size_t)ret >= sizeof(sim->name))
31 errx(1, "%s: could not write simulation name", __func__);
32 sim->ng = 0;
33 for (d = 0; d < 3; d++) {
34 sim->nd[d] = 1;
35 sim->origo[d] = 0.0;
36 sim->L[d] = 1.0;
37 }
38 sim->t = 0.0;
39 sim->t_end = 0.0;
40 sim->dt = 0.0;
41 sim->dt_file = 1.0;
42 sim->n_file = 0;
43 sim->iter = 0;
44 sim->t_file = 0;
45 sim->ng = 0;
46 sim->grains = NULL;
47 }
48
49 void
50 sim_free(struct simulation *sim)
51 {
52 free(sim->grains);
53 sim->grains = NULL;
54 sim->ng = 0;
55 }
56
57 void
58 sim_add_grain(struct simulation *sim, struct grain *g)
59 {
60 if (!(sim->grains = xreallocarray(sim->grains, sim->ng + 1, sizeof(*g))))
61 err(1, "%s: sim.grains reallocarray", __func__);
62 memcpy(&sim->grains[sim->ng++], g, sizeof(*g));
63 }
64
65 void
66 sim_read_grains(struct simulation *sim, FILE *stream)
67 {
68 char *line = NULL;
69 size_t linesize = 0;
70 ssize_t linelen;
71 struct grain *g;
72
73 while ((linelen = getline(&line, &linesize, stream)) > 0) {
74 g = grain_read(line);
75 sim_add_grain(sim, g);
76 free(g);
77 }
78
79 free(line);
80 }
81
82 void
83 sim_print_grains(const struct simulation *sim, FILE *stream)
84 {
85 grains_print(stream, sim->grains, sim->ng);
86 }
87
88 void
89 sim_print_grains_vtk(const struct simulation *sim, FILE *stream)
90 {
91 grains_print_vtk(stream, sim->grains, sim->ng);
92 }
93
94 void
95 sim_print_grains_energy(const struct simulation *sim, FILE *stream)
96 {
97 grains_print_energy(stream, sim->grains, sim->ng);
98 }
99
100 void
101 sim_write_output_files(struct simulation *sim)
102 {
103 int ret;
104 char outfile[256];
105 FILE *fp;
106
107 ret = snprintf(outfile, sizeof(outfile), "%s.grains.%05zu.tsv",
108 sim->name, sim->n_file++);
109 if (ret < 0 || (size_t)ret >= sizeof(outfile))
110 errx(1, "%s: could not write grains filename", __func__);
111 if ((fp = fopen(outfile, "w")) != NULL) {
112 sim_print_grains(sim, fp);
113 fclose(fp);
114 } else
115 err(1, "%s: fopen: %s", __func__, outfile);
116 }
117
118 double
119 sim_grain_pair_overlap(struct simulation *sim, size_t i, size_t j)
120 {
121 double pos[3], overlap;
122 int d;
123
124 for (d = 0; d < 3; d++)
125 pos[d] = sim->grains[i].pos[d] - sim->grains[j].pos[d];
126 overlap = (sim->grains[i].diameter + sim->grains[j].diameter) / 2.0
127 - euclidean_norm(pos, 3);
128
129 return overlap;
130 }
131
132 static void
133 sim_check_add_contact(struct simulation *sim, size_t i, size_t j)
134 {
135
136 double overlap, centerdist[3];
137 int d;
138
139 if (i >= j || !sim->grains[i].enabled || !sim->grains[j].enabled)
140 return;
141
142 for (d = 0; d < 3; d++)
143 centerdist[d] = sim->grains[i].pos[d] - sim->grains[j].pos[d];
144
145 overlap = 0.5*(sim->grains[i].diameter + sim->grains[j].diameter)
146 - euclidean_norm(centerdist, 3);
147
148 grain_register_contact(&sim->grains[i], i, j, centerdist, overlap);
149 }
150
151 void
152 sim_force_reset(struct simulation *sim)
153 {
154 size_t i;
155 for (i = 0; i < sim->ng; i++)
156 grain_force_reset(&sim->grains[i]);
157 }
158
159 void
160 sim_add_acceleration_scalar(struct simulation *sim, double acc, int dimension)
161 {
162 size_t i;
163
164 for (i = 0; i < sim->ng; i++)
165 sim->grains[i].forceext[dimension]
166 += acc * grain_mass(&sim->grains[i]);
167 }
168
169 /* Silbert et al. 2001 */
170 void
171 sim_set_timestep(struct simulation *sim)
172 {
173 int ret = 0;
174 size_t i;
175 double t_col, t_col_min = INFINITY;
176
177 for (i = 0; i < sim->ng; i++) {
178 t_col = grain_collision_time(&sim->grains[i]);
179 if (t_col_min > t_col)
180 t_col_min = t_col;
181 }
182
183 check_float("t_col_min", t_col_min, &ret);
184 sim->dt = t_col_min / 50.0;
185 }
186
187 /* TODO: add grid sorting */
188 void
189 sim_detect_contacts(struct simulation *sim)
190 {
191 size_t i, j;
192
193 for (i = 0; i < sim->ng; i++)
194 if (sim->grains[i].enabled)
195 for (j = 0; j < sim->ng; j++)
196 if (i < j)
197 sim_check_add_contact(sim, i, j);
198 }
199
200 void
201 sim_resolve_interactions(struct simulation *sim)
202 {
203 size_t i;
204 int ic;
205
206 for (i = 0; i < sim->ng; i++)
207 for (ic = 0; ic < MAXCONTACTS; ic++)
208 if (sim->grains[i].contacts[ic].active && i < sim->grains[i].contacts[ic].j)
209 grains_interact(&sim->grains[i],
210 &sim->grains[sim->grains[i].contacts[ic].j],
211 ic, sim->dt);
212 }
213
214 void
215 sim_step_time(struct simulation *sim)
216 {
217 size_t i;
218
219 sim_force_reset(sim);
220 sim_detect_contacts(sim);
221 sim_resolve_interactions(sim);
222
223 for (i = 0; i < sim->ng; i++)
224 grain_temporal_integration(&sim->grains[i], sim->dt);
225
226 if ((sim->t_file >= sim->dt_file || sim->iter == 0) &&
227 strncmp(sim->name, DEFAULT_SIMULATION_NAME,
228 sizeof(sim->name)) != 0) {
229 sim_write_output_files(sim);
230 sim->t_file = 0.0;
231 }
232
233 sim->t += sim->dt;
234 sim->t_file += sim->dt;
235 sim->iter++;
236 }
237
238 void
239 sim_run_time_loop(struct simulation *sim)
240 {
241 if (sim->dt > sim->dt_file)
242 sim->dt = sim->dt_file;
243
244 do {
245 sim_step_time(sim);
246 } while (sim->t - sim->dt < sim->t_end);
247 }