tsimulation.c - slidergrid - grid of elastic sliders on a frictional surface
HTML git clone git://src.adamsgaard.dk/slidergrid
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tsimulation.c (13479B)
---
1 #include <stdio.h>
2 #include <math.h>
3 #include "slider.h"
4 #include "simulation.h"
5 #include "constants.h"
6
7 simulation create_simulation()
8 {
9 simulation sim;
10 sim.N = 0;
11 sim.time = 0.0;
12 sim.time_end = 0.0;
13 sim.dt = 0.1;
14 sim.iteration = 0;
15 sim.bond_length_limit = 0;
16 sim.id = "unnamed";
17 sim.file_interval = 0.0;
18 sim.file_number = 0;
19 sim.verbose = 0;
20 return sim;
21 }
22
23 int check_simulation_values(const simulation* sim)
24 {
25 int return_status = 0;
26
27 // Slider-specific parameters
28 int i;
29 for (i=0; i<sim->N; i++) {
30
31 if (sim->sliders[i].mass <= 0.0) {
32 fprintf(stderr, "Error: Mass of slider %d is zero or negative "
33 "(%f kg)\n", i, sim->sliders[i].mass);
34 return_status = 1;
35 }
36
37 if (sim->sliders[i].moment_of_inertia <= 0.0) {
38 fprintf(stderr, "Error: Moment of inertia of slider %d is "
39 "zero or negative (%f kg*m*m)\n",
40 i, sim->sliders[i].moment_of_inertia);
41 return_status = 1;
42 }
43
44 }
45
46 // General parameters
47 if (sim->N <= 0) {
48 fprintf(stderr, "Error: The number of sliders (N = %d) must be a "
49 "positive number.\n", sim->N);
50 return_status = 1;
51 }
52
53 if (sim->dt <= 0.0) {
54 fprintf(stderr, "Error: The numerical time step (dt = %f) must be a "
55 "positive value.\n", sim->dt);
56 return_status = 1;
57 }
58
59 if (sim->dt > 1.0e9) {
60 fprintf(stderr, "Error: The numerical time step (dt = %f) is "
61 "invalid. Did you set bond_parallel_stiffness?\n", sim->dt);
62 return_status = 1;
63 }
64
65 if (sim->time > sim->time_end) {
66 fprintf(stderr, "Error: Current time (time = %f s) exceeds "
67 "total time (time_end = %f s)\n",
68 sim->time, sim->time_end);
69 return_status = 1;
70 }
71
72 if (sim->bond_length_limit <= 0.0) {
73 fprintf(stderr, "Error: The inter-slider bond length limit "
74 "(bond_length_limit = %f) must be a positive value.\n",
75 sim->bond_length_limit);
76 return_status = 1;
77 }
78
79 if (sim->file_interval < 0.0) {
80 fprintf(stderr, "Error: The output file interval"
81 "(file_interval = %f) must be a zero or positive value.\n",
82 sim->file_interval);
83 return_status = 1;
84 }
85
86 return return_status;
87 }
88
89 Float find_time_step(const slider* sliders, const int N, const Float safety)
90 {
91 Float smallest_mass = 1.0e20;
92 Float largest_stiffness = 0.0;
93
94 int i;
95 for (i=0; i<N; i++) {
96 if (sliders[i].mass < smallest_mass)
97 smallest_mass = sliders[i].mass;
98 if (sliders[i].bond_parallel_kv_stiffness > largest_stiffness)
99 largest_stiffness = sliders[i].bond_parallel_kv_stiffness;
100 }
101
102 return safety/sqrt(largest_stiffness/smallest_mass);
103 }
104
105 int save_slider_data_to_file(
106 const slider* sliders,
107 const int N,
108 const char* filename)
109 {
110 FILE* f = fopen(filename, "w");
111 if (f == NULL) {
112 fprintf(stderr, "Error: Could not open output file %s.", filename);
113 return 1;
114 }
115
116 int i;
117 for (i=0; i<N; i++)
118 fprintf(f,
119 "%f\t%f\t%f\t" // pos
120 "%f\t%f\t%f\t" // vel
121 "%f\t%f\t%f\t" // acc
122 "%f\t%f\t%f\t" // force
123 "%f\t%f\t%f\t" // angpos
124 "%f\t%f\t%f\t" // angvel
125 "%f\t%f\t%f\t" // angacc
126 "%f\t%f\t%f\t" // torque
127 "%f\t" // mass
128 "%f" // moment of inertia
129 "\n",
130 sliders[i].pos.x, sliders[i].pos.y, sliders[i].pos.z,
131 sliders[i].vel.x, sliders[i].vel.y, sliders[i].vel.z,
132 sliders[i].acc.x, sliders[i].acc.y, sliders[i].acc.z,
133 sliders[i].force.x, sliders[i].force.y, sliders[i].force.z,
134 sliders[i].angpos.x, sliders[i].angpos.y, sliders[i].angpos.z,
135 sliders[i].angvel.x, sliders[i].angvel.y, sliders[i].angvel.z,
136 sliders[i].angacc.x, sliders[i].angacc.y, sliders[i].angacc.z,
137 sliders[i].torque.x, sliders[i].torque.y, sliders[i].torque.z,
138 sliders[i].mass,
139 sliders[i].moment_of_inertia
140 );
141
142 fclose(f);
143 return 0;
144 }
145
146 int save_status_to_file(const simulation* sim, const char* filename)
147 {
148 FILE* f = fopen(filename, "w");
149 if (f == NULL) {
150 fprintf(stderr, "Error: Could not open output file %s.", filename);
151 return 1;
152 }
153
154 fprintf(f,
155 "%f\t" // current time
156 "%f\t" // percentage completed
157 "%d" // last output file number
158 ,
159 sim->time,
160 100.0*sim->time/sim->time_end,
161 sim->file_number
162 );
163
164 fclose(f);
165 return 0;
166
167 }
168
169 int save_general_state_to_file(const simulation* sim, const char* filename)
170 {
171 FILE* f = fopen(filename, "w");
172 if (f == NULL) {
173 fprintf(stderr, "Error: Could not open output file %s.", filename);
174 return 1;
175 }
176
177 fprintf(f,
178 "%f\t" // VERSION
179 "%s\t" // sim->id
180 "%d\t" // sim->N
181 "%f\t" // sim->time
182 "%f\t" // sim->time_end
183 "%f\t" // sim->dt
184 "%f\t" // sim->file_interval
185 "%ld\t" // sim->iteration
186 "%f" // sim->bond_length_limit
187 ,
188 VERSION,
189 sim->id,
190 sim->N,
191 sim->time,
192 sim->time_end,
193 sim->dt,
194 sim->file_interval,
195 sim->iteration,
196 sim->bond_length_limit
197 );
198
199 fclose(f);
200 return 0;
201 }
202
203 /* save all simulation information to binary format */
204 int save_binary_file(const simulation* sim, const char* filename)
205 {
206 FILE* f = fopen(filename, "wb");
207 if (f == NULL) {
208 fprintf(stderr, "Error: Could not open output file %s.", filename);
209 return 1;
210 }
211
212
213
214 fclose(f);
215 return 0;
216 }
217
218 /* save slider information as an unstructured grid to a VTK xml file. The
219 * filename extension should be '.vtu' in order for Paraview to correctly handle
220 * the file upon import. See VTK-file-formats.pdf for documentation on the file
221 * format. */
222 int save_sliders_to_vtk_file(
223 const slider* sliders,
224 const int N,
225 const char* filename)
226 {
227 FILE* f = fopen(filename, "w");
228 if (f == NULL) {
229 fprintf(stderr, "Error: Could not open output file %s.", filename);
230 return 1;
231 }
232 int i;
233
234 fprintf(f, "<?xml version=\"1.0\"?>\n");
235 fprintf(f, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
236 "byte_order=\"LittleEndian\">\n");
237 fprintf(f, " <UnstructuredGrid>\n");
238 fprintf(f, " <Piece NumberOfPoints=\"%d\" NumberOfCells=\"0\">\n", N);
239
240 fprintf(f, " <Points>\n");
241
242 // positions
243 fprintf(f, " <DataArray Name=\"Position [m]\" type=\"Float32\" "
244 "NumberOfComponents=\"3\" format=\"ascii\">\n");
245 fprintf(f, " ");
246 for (i=0; i<N; i++)
247 fprintf(f, "%f %f %f ",
248 sliders[i].pos.x, sliders[i].pos.y, sliders[i].pos.z);
249 fprintf(f, "\n </DataArray>\n");
250
251 fprintf(f, " </Points>\n");
252
253 fprintf(f, " <PointData Scalars=\"Mass [kg]\" "
254 "Vectors=\"vector\">\n");
255
256 // mass
257 fprintf(f, " <DataArray type=\"Float32\" "
258 "Name=\"Mass [kg]\" format=\"ascii\">\n");
259 fprintf(f, " ");
260 for (i=0; i<N; i++)
261 fprintf(f, "%f ", sliders[i].mass);
262 fprintf(f, "\n </DataArray>\n");
263
264 // velocities
265 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
266 "Name=\"Velocity [m/s]\" format=\"ascii\">\n");
267 fprintf(f, " ");
268 for (i=0; i<N; i++)
269 fprintf(f, "%f %f %f ",
270 sliders[i].vel.x, sliders[i].vel.y, sliders[i].vel.z);
271 fprintf(f, "\n </DataArray>\n");
272
273 // accelerations
274 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
275 "Name=\"Acceleration [m/(s*s)]\" format=\"ascii\">\n");
276 fprintf(f, " ");
277 for (i=0; i<N; i++)
278 fprintf(f, "%f %f %f ",
279 sliders[i].acc.x, sliders[i].acc.y, sliders[i].acc.z);
280 fprintf(f, "\n </DataArray>\n");
281
282 // force
283 fprintf(f, "\n <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
284 "Name=\"Force [N]\" format=\"ascii\">\n");
285 fprintf(f, " ");
286 for (i=0; i<N; i++)
287 fprintf(f, "%f %f %f ",
288 sliders[i].force.x, sliders[i].force.y, sliders[i].force.z);
289 fprintf(f, "\n </DataArray>\n");
290
291 // angular positions
292 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
293 "Name=\"Angular position [rad]\" format=\"ascii\">\n");
294 fprintf(f, " ");
295 for (i=0; i<N; i++)
296 fprintf(f, "%f %f %f ",
297 sliders[i].angpos.x, sliders[i].angpos.y, sliders[i].angpos.z);
298 fprintf(f, "\n </DataArray>\n");
299
300 // angular velocities
301 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
302 "Name=\"Angular velocity [rad/s]\" format=\"ascii\">\n");
303 fprintf(f, " ");
304 for (i=0; i<N; i++)
305 fprintf(f, "%f %f %f ",
306 sliders[i].angvel.x, sliders[i].angvel.y, sliders[i].angvel.z);
307 fprintf(f, "\n </DataArray>\n");
308
309 // angular accelerations
310 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
311 "Name=\"Angular acceleration [rad/(s*s)]\" format=\"ascii\">\n");
312 fprintf(f, " ");
313 for (i=0; i<N; i++)
314 fprintf(f, "%f %f %f ",
315 sliders[i].angacc.x, sliders[i].angacc.y, sliders[i].angacc.z);
316 fprintf(f, "\n </DataArray>\n");
317
318 // torque
319 fprintf(f, " <DataArray NumberOfComponents=\"3\" type=\"Float32\" "
320 "Name=\"Torque [N*m]\" format=\"ascii\">\n");
321 fprintf(f, " ");
322 for (i=0; i<N; i++)
323 fprintf(f, "%f %f %f ",
324 sliders[i].torque.x, sliders[i].torque.y, sliders[i].torque.z);
325 fprintf(f, "\n </DataArray>\n");
326
327 // moment of inertia
328 fprintf(f, " <DataArray type=\"Float32\" "
329 "Name=\"Moment of inertia [kg*m*m]\" format=\"ascii\">\n");
330 fprintf(f, " ");
331 for (i=0; i<N; i++)
332 fprintf(f, "%f ", sliders[i].moment_of_inertia);
333 fprintf(f, "\n </DataArray>\n");
334
335 // bond-parallel stiffness
336 fprintf(f, " <DataArray type=\"Float32\" "
337 "Name=\"Bond-parallel stiffness [N/m]\" format=\"ascii\">\n");
338 fprintf(f, "\n ");
339 for (i=0; i<N; i++)
340 fprintf(f, "%f ", sliders[i].bond_parallel_kv_stiffness);
341 fprintf(f, " </DataArray>\n");
342
343 // bond-parallel viscosity
344 fprintf(f, " <DataArray type=\"Float32\" "
345 "Name=\"Bond-parallel viscosity [N/(m/s)]\" format=\"ascii\">\n");
346 fprintf(f, " ");
347 for (i=0; i<N; i++)
348 fprintf(f, "%f ", sliders[i].bond_parallel_kv_viscosity);
349 fprintf(f, "\n </DataArray>\n");
350
351 fprintf(f, " </PointData>\n");
352 fprintf(f, " <CellData>\n");
353 fprintf(f, " </CellData>\n");
354 fprintf(f, " <Cells>\n");
355 fprintf(f, " <DataArray type=\"Int32\" Name=\"connectivity\" "
356 "format=\"ascii\">\n");
357 fprintf(f, " </DataArray>\n");
358 fprintf(f, " <DataArray type=\"Int32\" Name=\"offsets\" "
359 "format=\"ascii\">\n");
360 fprintf(f, " </DataArray>\n");
361 fprintf(f, " <DataArray type=\"UInt8\" Name=\"types\" "
362 "format=\"ascii\">\n");
363 fprintf(f, " </DataArray>\n");
364 fprintf(f, " </Cells>\n");
365 fprintf(f, " </Piece>\n");
366 fprintf(f, " </UnstructuredGrid>\n");
367 fprintf(f, "</VTKFile>\n");
368
369 fclose(f);
370 return 0;
371 }
372
373 int write_simulation_output(simulation* sim, char* output_folder)
374 {
375 char filename[1000];
376
377 // slider parameters
378 sprintf(filename, "%s/%s.sliders.%06d.txt",
379 output_folder, sim->id, sim->file_number);
380 if (save_slider_data_to_file(sim->sliders, sim->N, filename)) {
381 fprintf(stderr, "\nFatal error: Could not save to output file "
382 "'%s'.\n", filename);
383 return 1;
384 }
385
386 // other parameters
387 sprintf(filename, "%s/%s.general.%06d.txt",
388 output_folder, sim->id, sim->file_number);
389 if (save_general_state_to_file(sim, filename)) {
390 fprintf(stderr, "\nFatal error: Could not save to output file "
391 "'%s'.\n", filename);
392 return 1;
393 }
394
395 // sliders to VTK file
396 sprintf(filename, "%s/%s.sliders.%06d.vtu",
397 output_folder, sim->id, sim->file_number);
398 if (save_sliders_to_vtk_file(sim->sliders, sim->N, filename)) {
399 fprintf(stderr, "\nFatal error: Could not save to output file "
400 "'%s'.\n", filename);
401 return 1;
402 }
403
404 // simulation status, temporal values and last output file number
405 sprintf(filename, "%s/%s.status.dat", output_folder, sim->id);
406 if (save_status_to_file(sim, filename)) {
407 fprintf(stderr, "\nFatal error: Could not save to output file "
408 "'%s'.\n", filename);
409 return 1;
410 }
411
412 sim->file_number++;
413 return 0;
414 }