core.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
core.py (34926B)
---
1 import numpy
2 from .common import VERSION
3 from .analysis import SimAnalysis
4 from .fileio import SimIO
5 from .fluid import SimFluid
6 from .plotting import SimPlotting
7 from .runner import SimRun
8 from .visualize import SimVisualize
9 from .world import SimSetup
10
11
12 class sim(SimIO, SimSetup, SimFluid, SimAnalysis, SimPlotting,
13 SimVisualize, SimRun):
14 '''
15 Class containing all ``sphere`` data.
16
17 Contains functions for reading and writing binaries, as well as simulation
18 setup and data analysis. Most arrays are initialized to default values.
19
20 :param np: The number of particles to allocate memory for (default=1)
21 :type np: int
22 :param nd: The number of spatial dimensions (default=3). Note that 2D and
23 1D simulations currently are not possible.
24 :type nd: int
25 :param nw: The number of dynamic walls (default=1)
26 :type nw: int
27 :param sid: The simulation id (default='unnamed'). The simulation files
28 will be written with this base name.
29 :type sid: str
30 :param fluid: Setup fluid simulation (default=False)
31 :type fluid: bool
32 :param cfd_solver: Fluid solver to use if fluid == True. 0: Navier-Stokes
33 (default), 1: Darcy.
34 :type cfd_solver: int
35 '''
36
37 def __init__(self, sid='unnamed', np=0, nd=3, nw=0, fluid=False):
38
39 # Sphere version number
40 self.version = numpy.ones(1, dtype=numpy.float64)*VERSION
41
42 # The number of spatial dimensions. Values other that 3 do not work
43 self.nd = int(nd)
44
45 # The number of particles
46 self.np = int(np)
47
48 # The simulation id (text string)
49 self.sid = sid
50
51 ## Time parameters
52 # Computational time step length [s]
53 self.time_dt = numpy.zeros(1, dtype=numpy.float64)
54
55 # Current time [s]
56 self.time_current = numpy.zeros(1, dtype=numpy.float64)
57
58 # Total time [s]
59 self.time_total = numpy.zeros(1, dtype=numpy.float64)
60
61 # File output interval [s]
62 self.time_file_dt = numpy.zeros(1, dtype=numpy.float64)
63
64 # The number of files written
65 self.time_step_count = numpy.zeros(1, dtype=numpy.uint32)
66
67 ## World dimensions and grid data
68 # The Euclidean coordinate to the origo of the sorting grid
69 self.origo = numpy.zeros(self.nd, dtype=numpy.float64)
70
71 # The sorting grid size (x, y, z)
72 self.L = numpy.zeros(self.nd, dtype=numpy.float64)
73
74 # The number of sorting cells in each dimension
75 self.num = numpy.zeros(self.nd, dtype=numpy.uint32)
76
77 # Whether to treat the lateral boundaries as periodic (1) or not (0)
78 self.periodic = numpy.zeros(1, dtype=numpy.uint32)
79
80 # Adaptively resize grid to assemblage height (0: no, 1: yes)
81 self.adaptive = numpy.zeros(1, dtype=numpy.uint32)
82
83 ## Particle data
84 # Particle position vectors [m]
85 self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
86
87 # Particle radii [m]
88 self.radius = numpy.ones(self.np, dtype=numpy.float64)
89
90 # The sums of x and y movement [m]
91 self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64)
92
93 # The linear velocities [m/s]
94 self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
95
96 # Fix the particle kinematics?
97 # 0: No (DEFAULT, don't fix linear or angular acceleration)
98 # 1: Yes (fix horizontal movement, allow vertical movement, disable rotation)
99 # 10: Yes (fix horizontal movement, allow vertical movement, disable rotation)
100 # -1: Yes (fix all linear and rotational movement)
101 # -10: Yes (fix all rotational movement)
102 self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
103
104 # The linear force vectors [N]
105 self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
106
107 # The angular position vectors [rad]
108 self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
109
110 # The angular velocity vectors [rad/s]
111 self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
112
113 # The torque vectors [N*m]
114 self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
115
116 # The shear friction energy dissipation rates [W]
117 self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
118
119 # The total shear energy dissipations [J]
120 self.es = numpy.zeros(self.np, dtype=numpy.float64)
121
122 # The viscous energy dissipation rates [W]
123 self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
124
125 # The total viscois energy dissipation [J]
126 self.ev = numpy.zeros(self.np, dtype=numpy.float64)
127
128 # The total particle pressures [Pa]
129 self.p = numpy.zeros(self.np, dtype=numpy.float64)
130
131 # The gravitational acceleration vector [N*m/s]
132 self.g = numpy.array([0.0, 0.0, 0.0], dtype=numpy.float64)
133
134 # The Hookean coefficient for elastic stiffness normal to the contacts
135 # [N/m]
136 self.k_n = numpy.ones(1, dtype=numpy.float64) * 1.16e9
137
138 # The Hookean coefficient for elastic stiffness tangential to the
139 # contacts [N/m]
140 self.k_t = numpy.ones(1, dtype=numpy.float64) * 1.16e9
141
142 # The Hookean coefficient for elastic stiffness opposite of contact
143 # rotations. UNUSED
144 self.k_r = numpy.zeros(1, dtype=numpy.float64)
145
146 # Young's modulus for contact stiffness [Pa]. This value is used
147 # instead of the Hookean stiffnesses (k_n, k_t) when self.E is larger
148 # than 0.0.
149 self.E = numpy.zeros(1, dtype=numpy.float64)
150
151 # The viscosity normal to the contact [N/(m/s)]
152 self.gamma_n = numpy.zeros(1, dtype=numpy.float64)
153
154 # The viscosity tangential to the contact [N/(m/s)]
155 self.gamma_t = numpy.zeros(1, dtype=numpy.float64)
156
157 # The viscosity to contact rotation [N/(m/s)]
158 self.gamma_r = numpy.zeros(1, dtype=numpy.float64)
159
160 # The coefficient of static friction on the contact [-]
161 self.mu_s = numpy.ones(1, dtype=numpy.float64) * 0.5
162
163 # The coefficient of dynamic friction on the contact [-]
164 self.mu_d = numpy.ones(1, dtype=numpy.float64) * 0.5
165
166 # The coefficient of rotational friction on the contact [-]
167 self.mu_r = numpy.zeros(1, dtype=numpy.float64)
168
169 # The viscosity normal to the walls [N/(m/s)]
170 self.gamma_wn = numpy.zeros(1, dtype=numpy.float64)
171
172 # The viscosity tangential to the walls [N/(m/s)]
173 self.gamma_wt = numpy.zeros(1, dtype=numpy.float64)
174
175 # The coeffient of static friction of the walls [-]
176 self.mu_ws = numpy.ones(1, dtype=numpy.float64) * 0.5
177
178 # The coeffient of dynamic friction of the walls [-]
179 self.mu_wd = numpy.ones(1, dtype=numpy.float64) * 0.5
180
181 # The particle density [kg/(m^3)]
182 self.rho = numpy.ones(1, dtype=numpy.float64) * 2600.0
183
184 # The contact model to use
185 # 1: Normal: elasto-viscous, tangential: visco-frictional
186 # 2: Normal: elasto-viscous, tangential: elasto-visco-frictional
187 self.contactmodel = numpy.ones(1, dtype=numpy.uint32) * 2 # lin-visc-el
188
189 # Capillary bond prefactor
190 self.kappa = numpy.zeros(1, dtype=numpy.float64)
191
192 # Capillary bond debonding distance [m]
193 self.db = numpy.zeros(1, dtype=numpy.float64)
194
195 # Capillary bond liquid volume [m^3]
196 self.V_b = numpy.zeros(1, dtype=numpy.float64)
197
198 ## Wall data
199 # Number of dynamic walls
200 # nw=1: Uniaxial (also used for shear experiments)
201 # nw=2: Biaxial
202 # nw=5: Triaxial
203 self.nw = int(nw)
204
205 # Wall modes
206 # 0: Fixed
207 # 1: Normal stress condition
208 # 2: Normal velocity condition
209 # 3: Normal stress and shear stress condition
210 self.wmode = numpy.zeros(self.nw, dtype=numpy.int32)
211
212 # Wall normals
213 self.w_n = numpy.zeros((self.nw, self.nd), dtype=numpy.float64)
214 if self.nw >= 1:
215 self.w_n[0, 2] = -1.0
216 if self.nw >= 2:
217 self.w_n[1, 0] = -1.0
218 if self.nw >= 3:
219 self.w_n[2, 0] = 1.0
220 if self.nw >= 4:
221 self.w_n[3, 1] = -1.0
222 if self.nw >= 5:
223 self.w_n[4, 1] = 1.0
224
225 # Wall positions on the axes that are parallel to the wall normal [m]
226 self.w_x = numpy.ones(self.nw, dtype=numpy.float64)
227
228 # Wall masses [kg]
229 self.w_m = numpy.zeros(self.nw, dtype=numpy.float64)
230
231 # Wall velocities on the axes that are parallel to the wall normal [m/s]
232 self.w_vel = numpy.zeros(self.nw, dtype=numpy.float64)
233
234 # Wall forces on the axes that are parallel to the wall normal [m/s]
235 self.w_force = numpy.zeros(self.nw, dtype=numpy.float64)
236
237 # Wall stress on the axes that are parallel to the wall normal [Pa]
238 self.w_sigma0 = numpy.zeros(self.nw, dtype=numpy.float64)
239
240 # Wall stress modulation amplitude [Pa]
241 self.w_sigma0_A = numpy.zeros(1, dtype=numpy.float64)
242
243 # Wall stress modulation frequency [Hz]
244 self.w_sigma0_f = numpy.zeros(1, dtype=numpy.float64)
245
246 # Wall shear stress, enforced when wmode == 3
247 self.w_tau_x = numpy.zeros(1, dtype=numpy.float64)
248
249 ## Bond parameters
250 # Radius multiplier to the parallel-bond radii
251 self.lambda_bar = numpy.ones(1, dtype=numpy.float64)
252
253 # Number of bonds
254 self.nb0 = 0
255
256 # Bond tensile strength [Pa]
257 self.sigma_b = numpy.ones(1, dtype=numpy.float64) * numpy.inf
258
259 # Bond shear strength [Pa]
260 self.tau_b = numpy.ones(1, dtype=numpy.float64) * numpy.inf
261
262 # Bond pairs
263 self.bonds = numpy.zeros((self.nb0, 2), dtype=numpy.uint32)
264
265 # Parallel bond movement
266 self.bonds_delta_n = numpy.zeros(self.nb0, dtype=numpy.float64)
267
268 # Shear bond movement
269 self.bonds_delta_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
270
271 # Twisting bond movement
272 self.bonds_omega_n = numpy.zeros(self.nb0, dtype=numpy.float64)
273
274 # Bending bond movement
275 self.bonds_omega_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
276
277 ## Fluid parameters
278
279 # Simulate fluid? True: Yes, False: no
280 self.fluid = fluid
281
282 if self.fluid:
283
284 # Fluid solver type
285 # 0: Navier Stokes (fluid with inertia)
286 # 1: Stokes-Darcy (fluid without inertia)
287 self.cfd_solver = numpy.zeros(1, dtype=numpy.int32)
288
289 # Fluid dynamic viscosity [N/(m/s)]
290 self.mu = numpy.zeros(1, dtype=numpy.float64)
291
292 # Fluid velocities [m/s]
293 self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
294 dtype=numpy.float64)
295
296 # Fluid pressures [Pa]
297 self.p_f = numpy.zeros((self.num[0], self.num[1], self.num[2]),
298 dtype=numpy.float64)
299
300 # Fluid cell porosities [-]
301 self.phi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
302 dtype=numpy.float64)
303
304 # Fluid cell porosity change [1/s]
305 self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
306 dtype=numpy.float64)
307
308 # Fluid density [kg/(m^3)]
309 self.rho_f = numpy.ones(1, dtype=numpy.float64) * 1.0e3
310
311 # Pressure modulation at the top boundary
312 self.p_mod_A = numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
313 self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
314 self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
315
316 ## Fluid solver parameters
317
318 if self.cfd_solver[0] == 1: # Darcy solver
319 # Boundary conditions at the sides of the fluid grid
320 # 0: Dirichlet
321 # 1: Neumann
322 # 2: Periodic (default)
323 self.bc_xn = numpy.ones(1, dtype=numpy.int32)*2 # Neg. x bc
324 self.bc_xp = numpy.ones(1, dtype=numpy.int32)*2 # Pos. x bc
325 self.bc_yn = numpy.ones(1, dtype=numpy.int32)*2 # Neg. y bc
326 self.bc_yp = numpy.ones(1, dtype=numpy.int32)*2 # Pos. y bc
327
328 # Boundary conditions at the top and bottom of the fluid grid
329 # 0: Dirichlet (default)
330 # 1: Neumann free slip
331 # 2: Neumann no slip (Navier Stokes), Periodic (Darcy)
332 # 3: Periodic (Navier-Stokes solver only)
333 # 4: Constant flux (Darcy solver only)
334 self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
335 self.bc_top = numpy.zeros(1, dtype=numpy.int32)
336 # Free slip boundaries? 1: yes
337 self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
338 self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
339
340 # Boundary-normal flux (in case of bc_*=4)
341 self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
342 self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
343
344 # Hold pressures constant in fluid cell (0: True, 1: False)
345 self.p_f_constant = numpy.zeros((self.num[0],
346 self.num[1],
347 self.num[2]), dtype=numpy.int32)
348
349 # Navier-Stokes
350 if self.cfd_solver[0] == 0:
351
352 # Smoothing parameter, should be in the range [0.0;1.0[.
353 # 0.0=no smoothing.
354 self.gamma = numpy.array(0.0)
355
356 # Under-relaxation parameter, should be in the range ]0.0;1.0].
357 # 1.0=no under-relaxation
358 self.theta = numpy.array(1.0)
359
360 # Velocity projection parameter, should be in the range
361 # [0.0;1.0]
362 self.beta = numpy.array(0.0)
363
364 # Tolerance criteria for the normalized max. residual
365 self.tolerance = numpy.array(1.0e-3)
366
367 # The maximum number of iterations to perform per time step
368 self.maxiter = numpy.array(1e4)
369
370 # The number of DEM time steps to perform between CFD updates
371 self.ndem = numpy.array(1)
372
373 # Porosity scaling factor
374 self.c_phi = numpy.ones(1, dtype=numpy.float64)
375
376 # Fluid velocity scaling factor
377 self.c_v = numpy.ones(1, dtype=numpy.float64)
378
379 # DEM-CFD time scaling factor
380 self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
381
382 ## Interaction forces
383 self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
384 self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
385 self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
386 self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
387
388 # Darcy
389 elif self.cfd_solver[0] == 1:
390
391 # Tolerance criteria for the normalized max. residual
392 self.tolerance = numpy.array(1.0e-3)
393
394 # The maximum number of iterations to perform per time step
395 self.maxiter = numpy.array(1e4)
396
397 # The number of DEM time steps to perform between CFD updates
398 self.ndem = numpy.array(1)
399
400 # Porosity scaling factor
401 self.c_phi = numpy.ones(1, dtype=numpy.float64)
402
403 # Interaction forces
404 self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
405
406 # Adiabatic fluid compressibility [1/Pa].
407 # Fluid bulk modulus=1/self.beta_f
408 self.beta_f = numpy.ones(1, dtype=numpy.float64)*4.5e-10
409
410 # Hydraulic permeability prefactor [m*m]
411 self.k_c = numpy.ones(1, dtype=numpy.float64)*4.6e-10
412
413 else:
414 raise Exception('Value of cfd_solver not understood (' + \
415 str(self.cfd_solver[0]) + ')')
416
417 # Particle color marker
418 self.color = numpy.zeros(self.np, dtype=numpy.int32)
419
420 def __eq__(self, other):
421 '''
422 Called when to sim objects are compared. Returns 0 if the values
423 are identical.
424 '''
425 if self.version != other.version:
426 print('version')
427 return False
428 elif self.nd != other.nd:
429 print('nd')
430 return False
431 elif self.np != other.np:
432 print('np')
433 return False
434 elif self.time_dt != other.time_dt:
435 print('time_dt')
436 return False
437 elif self.time_current != other.time_current:
438 print('time_current')
439 return False
440 elif self.time_total != other.time_total:
441 print('time_total')
442 return False
443 elif self.time_file_dt != other.time_file_dt:
444 print('time_file_dt')
445 return False
446 elif self.time_step_count != other.time_step_count:
447 print('time_step_count')
448 return False
449 elif (self.origo != other.origo).any():
450 print('origo')
451 return False
452 elif (self.L != other.L).any():
453 print('L')
454 return 11
455 elif (self.num != other.num).any():
456 print('num')
457 return False
458 elif self.periodic != other.periodic:
459 print('periodic')
460 return False
461 elif self.adaptive != other.adaptive:
462 print('adaptive')
463 return False
464 elif (self.x != other.x).any():
465 print('x')
466 return False
467 elif (self.radius != other.radius).any():
468 print('radius')
469 return False
470 elif (self.xyzsum != other.xyzsum).any():
471 print('xyzsum')
472 return False
473 elif (self.vel != other.vel).any():
474 print('vel')
475 return False
476 elif (self.fixvel != other.fixvel).any():
477 print('fixvel')
478 return False
479 elif (self.force != other.force).any():
480 print('force')
481 return False
482 elif (self.angpos != other.angpos).any():
483 print('angpos')
484 return False
485 elif (self.angvel != other.angvel).any():
486 print('angvel')
487 return False
488 elif (self.torque != other.torque).any():
489 print('torque')
490 return False
491 elif (self.es_dot != other.es_dot).any():
492 print('es_dot')
493 return False
494 elif (self.es != other.es).any():
495 print('es')
496 return False
497 elif (self.ev_dot != other.ev_dot).any():
498 print('ev_dot')
499 return False
500 elif (self.ev != other.ev).any():
501 print('ev')
502 return False
503 elif (self.p != other.p).any():
504 print('p')
505 return False
506 elif (self.g != other.g).any():
507 print('g')
508 return False
509 elif self.k_n != other.k_n:
510 print('k_n')
511 return False
512 elif self.k_t != other.k_t:
513 print('k_t')
514 return False
515 elif self.k_r != other.k_r:
516 print('k_r')
517 return False
518 elif self.E != other.E:
519 print('E')
520 return False
521 elif self.gamma_n != other.gamma_n:
522 print('gamma_n')
523 return False
524 elif self.gamma_t != other.gamma_t:
525 print('gamma_t')
526 return False
527 elif self.gamma_r != other.gamma_r:
528 print('gamma_r')
529 return False
530 elif self.mu_s != other.mu_s:
531 print('mu_s')
532 return False
533 elif self.mu_d != other.mu_d:
534 print('mu_d')
535 return False
536 elif self.mu_r != other.mu_r:
537 print('mu_r')
538 return False
539 elif self.rho != other.rho:
540 print('rho')
541 return False
542 elif self.contactmodel != other.contactmodel:
543 print('contactmodel')
544 return False
545 elif self.kappa != other.kappa:
546 print('kappa')
547 return False
548 elif self.db != other.db:
549 print('db')
550 return False
551 elif self.V_b != other.V_b:
552 print('V_b')
553 return False
554 elif self.nw != other.nw:
555 print('nw')
556 return False
557 elif (self.wmode != other.wmode).any():
558 print('wmode')
559 return False
560 elif (self.w_n != other.w_n).any():
561 print('w_n')
562 return False
563 elif (self.w_x != other.w_x).any():
564 print('w_x')
565 return False
566 elif (self.w_m != other.w_m).any():
567 print('w_m')
568 return False
569 elif (self.w_vel != other.w_vel).any():
570 print('w_vel')
571 return False
572 elif (self.w_force != other.w_force).any():
573 print('w_force')
574 return False
575 elif (self.w_sigma0 != other.w_sigma0).any():
576 print('w_sigma0')
577 return False
578 elif self.w_sigma0_A != other.w_sigma0_A:
579 print('w_sigma0_A')
580 return False
581 elif self.w_sigma0_f != other.w_sigma0_f:
582 print('w_sigma0_f')
583 return False
584 elif self.w_tau_x != other.w_tau_x:
585 print('w_tau_x')
586 return False
587 elif self.gamma_wn != other.gamma_wn:
588 print('gamma_wn')
589 return False
590 elif self.gamma_wt != other.gamma_wt:
591 print('gamma_wt')
592 return False
593 elif self.lambda_bar != other.lambda_bar:
594 print('lambda_bar')
595 return False
596 elif self.nb0 != other.nb0:
597 print('nb0')
598 return False
599 elif self.sigma_b != other.sigma_b:
600 print('sigma_b')
601 return False
602 elif self.tau_b != other.tau_b:
603 print('tau_b')
604 return False
605 elif (self.bonds != other.bonds).any():
606 print('bonds')
607 return False
608 elif (self.bonds_delta_n != other.bonds_delta_n).any():
609 print('bonds_delta_n')
610 return False
611 elif (self.bonds_delta_t != other.bonds_delta_t).any():
612 print('bonds_delta_t')
613 return False
614 elif (self.bonds_omega_n != other.bonds_omega_n).any():
615 print('bonds_omega_n')
616 return False
617 elif (self.bonds_omega_t != other.bonds_omega_t).any():
618 print('bonds_omega_t')
619 return False
620 elif self.fluid != other.fluid:
621 print('fluid')
622 return False
623
624 if self.fluid:
625 if self.cfd_solver != other.cfd_solver:
626 print('cfd_solver')
627 return False
628 elif self.mu != other.mu:
629 print('mu')
630 return False
631 elif (self.v_f != other.v_f).any():
632 print('v_f')
633 return False
634 elif (self.p_f != other.p_f).any():
635 print('p_f')
636 return False
637 #elif self.phi != other.phi).any():
638 #return False # Porosities not initialized correctly
639 elif (self.dphi != other.dphi).any():
640 print('d_phi')
641 return False
642 elif self.rho_f != other.rho_f:
643 print('rho_f')
644 return False
645 elif self.p_mod_A != other.p_mod_A:
646 print('p_mod_A')
647 return False
648 elif self.p_mod_f != other.p_mod_f:
649 print('p_mod_f')
650 return False
651 elif self.p_mod_phi != other.p_mod_phi:
652 print('p_mod_phi')
653 return False
654 elif self.bc_bot != other.bc_bot:
655 print('bc_bot')
656 return False
657 elif self.bc_top != other.bc_top:
658 print('bc_top')
659 return False
660 elif self.free_slip_bot != other.free_slip_bot:
661 print('free_slip_bot')
662 return False
663 elif self.free_slip_top != other.free_slip_top:
664 print('free_slip_top')
665 return False
666 elif self.bc_bot_flux != other.bc_bot_flux:
667 print('bc_bot_flux')
668 return False
669 elif self.bc_top_flux != other.bc_top_flux:
670 print('bc_top_flux')
671 return False
672 elif (self.p_f_constant != other.p_f_constant).any():
673 print('p_f_constant')
674 return False
675
676 if self.cfd_solver == 0:
677 if self.gamma != other.gamma:
678 print('gamma')
679 return False
680 elif self.theta != other.theta:
681 print('theta')
682 return False
683 elif self.beta != other.beta:
684 print('beta')
685 return False
686 elif self.tolerance != other.tolerance:
687 print('tolerance')
688 return False
689 elif self.maxiter != other.maxiter:
690 print('maxiter')
691 return False
692 elif self.ndem != other.ndem:
693 print('ndem')
694 return False
695 elif self.c_phi != other.c_phi:
696 print('c_phi')
697 return 84
698 elif self.c_v != other.c_v:
699 print('c_v')
700 elif self.dt_dem_fac != other.dt_dem_fac:
701 print('dt_dem_fac')
702 return 85
703 elif (self.f_d != other.f_d).any():
704 print('f_d')
705 return 86
706 elif (self.f_p != other.f_p).any():
707 print('f_p')
708 return 87
709 elif (self.f_v != other.f_v).any():
710 print('f_v')
711 return 88
712 elif (self.f_sum != other.f_sum).any():
713 print('f_sum')
714 return 89
715
716 if self.cfd_solver == 1:
717 if self.tolerance != other.tolerance:
718 print('tolerance')
719 return False
720 elif self.maxiter != other.maxiter:
721 print('maxiter')
722 return False
723 elif self.ndem != other.ndem:
724 print('ndem')
725 return False
726 elif self.c_phi != other.c_phi:
727 print('c_phi')
728 return 84
729 elif (self.f_p != other.f_p).any():
730 print('f_p')
731 return 86
732 elif self.beta_f != other.beta_f:
733 print('beta_f')
734 return 87
735 elif self.k_c != other.k_c:
736 print('k_c')
737 return 88
738 elif self.bc_xn != other.bc_xn:
739 print('bc_xn')
740 return False
741 elif self.bc_xp != other.bc_xp:
742 print('bc_xp')
743 return False
744 elif self.bc_yn != other.bc_yn:
745 print('bc_yn')
746 return False
747 elif self.bc_yp != other.bc_yp:
748 print('bc_yp')
749 return False
750
751 if (self.color != other.color).any():
752 print('color')
753 return False
754
755 # All equal
756 return True
757
758 def id(self, sid=''):
759 '''
760 Returns or sets the simulation id/name, which is used to identify
761 simulation files in the output folders.
762
763 :param sid: The desired simulation id. If left blank the current
764 simulation id will be returned.
765 :type sid: str
766 :returns: The current simulation id if no new value is set.
767 :return type: str
768 '''
769 if sid == '':
770 return self.sid
771 else:
772 self.sid = sid
773
774 def idAppend(self, string):
775 '''
776 Append a string to the simulation id/name, which is used to identify
777 simulation files in the output folders.
778
779 :param string: The string to append to the simulation id (`self.sid`).
780 :type string: str
781 '''
782 self.sid += string
783
784 def addParticle(self, x, radius, xyzsum=numpy.zeros(3), vel=numpy.zeros(3),
785 fixvel=numpy.zeros(1), force=numpy.zeros(3),
786 angpos=numpy.zeros(3), angvel=numpy.zeros(3),
787 torque=numpy.zeros(3), es_dot=numpy.zeros(1),
788 es=numpy.zeros(1), ev_dot=numpy.zeros(1),
789 ev=numpy.zeros(1), p=numpy.zeros(1), color=0):
790 '''
791 Add a single particle to the simulation object. The only required
792 parameters are the position (x) and the radius (radius).
793
794 :param x: A vector pointing to the particle center coordinate.
795 :type x: numpy.array
796 :param radius: The particle radius
797 :type radius: float
798 :param vel: The particle linear velocity (default=[0, 0, 0])
799 :type vel: numpy.array
800 :param fixvel: 0: Do not fix particle velocity (default), 1: Fix
801 horizontal linear velocity, -1: Fix horizontal and vertical linear
802 velocity
803 :type fixvel: float
804 :param angpos: The particle angular position (default=[0, 0, 0])
805 :type angpos: numpy.array
806 :param angvel: The particle angular velocity (default=[0, 0, 0])
807 :type angvel: numpy.array
808 :param torque: The particle torque (default=[0, 0, 0])
809 :type torque: numpy.array
810 :param es_dot: The particle shear energy loss rate (default=0)
811 :type es_dot: float
812 :param es: The particle shear energy loss (default=0)
813 :type es: float
814 :param ev_dot: The particle viscous energy rate loss (default=0)
815 :type ev_dot: float
816 :param ev: The particle viscous energy loss (default=0)
817 :type ev: float
818 :param p: The particle pressure (default=0)
819 :type p: float
820 '''
821
822 self.np += 1
823
824 self.x = numpy.append(self.x, [x], axis=0)
825 self.radius = numpy.append(self.radius, radius)
826 self.vel = numpy.append(self.vel, [vel], axis=0)
827 self.xyzsum = numpy.append(self.xyzsum, [xyzsum], axis=0)
828 self.fixvel = numpy.append(self.fixvel, fixvel)
829 self.force = numpy.append(self.force, [force], axis=0)
830 self.angpos = numpy.append(self.angpos, [angpos], axis=0)
831 self.angvel = numpy.append(self.angvel, [angvel], axis=0)
832 self.torque = numpy.append(self.torque, [torque], axis=0)
833 self.es_dot = numpy.append(self.es_dot, es_dot)
834 self.es = numpy.append(self.es, es)
835 self.ev_dot = numpy.append(self.ev_dot, ev_dot)
836 self.ev = numpy.append(self.ev, ev)
837 self.p = numpy.append(self.p, p)
838 self.color = numpy.append(self.color, color)
839 if self.fluid:
840 self.f_d = numpy.append(self.f_d, [numpy.zeros(3)], axis=0)
841 self.f_p = numpy.append(self.f_p, [numpy.zeros(3)], axis=0)
842 self.f_v = numpy.append(self.f_v, [numpy.zeros(3)], axis=0)
843 self.f_sum = numpy.append(self.f_sum, [numpy.zeros(3)], axis=0)
844
845 def deleteParticle(self, i):
846 '''
847 Delete particle(s) with index ``i``.
848
849 :param i: One or more particle indexes to delete
850 :type i: int, list or numpy.array
851 '''
852
853 # The user wants to delete several particles, indexes in a numpy.array
854 if type(i) == numpy.ndarray:
855 self.np -= i.size
856
857 # The user wants to delete several particles, indexes in a Python list
858 elif type(i) == list:
859 self.np -= len(i)
860
861 # The user wants to delete a single particle with a integer index
862 else:
863 self.np -= 1
864
865 if type(i) == tuple:
866 raise Exception('Cannot parse tuples as index value. ' +
867 'Valid types are int, list and numpy.ndarray')
868
869
870 self.x = numpy.delete(self.x, i, axis=0)
871 self.radius = numpy.delete(self.radius, i)
872 self.vel = numpy.delete(self.vel, i, axis=0)
873 self.xyzsum = numpy.delete(self.xyzsum, i, axis=0)
874 self.fixvel = numpy.delete(self.fixvel, i)
875 self.force = numpy.delete(self.force, i, axis=0)
876 self.angpos = numpy.delete(self.angpos, i, axis=0)
877 self.angvel = numpy.delete(self.angvel, i, axis=0)
878 self.torque = numpy.delete(self.torque, i, axis=0)
879 self.es_dot = numpy.delete(self.es_dot, i)
880 self.es = numpy.delete(self.es, i)
881 self.ev_dot = numpy.delete(self.ev_dot, i)
882 self.ev = numpy.delete(self.ev, i)
883 self.p = numpy.delete(self.p, i)
884 self.color = numpy.delete(self.color, i)
885 if self.fluid:
886 # Darcy and Navier-Stokes
887 self.f_p = numpy.delete(self.f_p, i, axis=0)
888 if self.cfd_solver[0] == 0: # Navier-Stokes
889 self.f_d = numpy.delete(self.f_d, i, axis=0)
890 self.f_v = numpy.delete(self.f_v, i, axis=0)
891 self.f_sum = numpy.delete(self.f_sum, i, axis=0)
892
893 def deleteAllParticles(self):
894 '''
895 Deletes all particles in the simulation object.
896 '''
897 self.np = 0
898 self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
899 self.radius = numpy.ones(self.np, dtype=numpy.float64)
900 self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64)
901 self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
902 self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
903 self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
904 self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
905 self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
906 self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
907 self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
908 self.es = numpy.zeros(self.np, dtype=numpy.float64)
909 self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
910 self.ev = numpy.zeros(self.np, dtype=numpy.float64)
911 self.p = numpy.zeros(self.np, dtype=numpy.float64)
912 self.color = numpy.zeros(self.np, dtype=numpy.int32)
913 self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
914 self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
915 self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
916 self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)