tgetting_started.md - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tgetting_started.md (18044B)
---
1 # Getting started
2 In this section, it is assumed that [Julia](https://julialang.org) and
3 [Granular.jl](https://github.com/anders-dc/Granular.jl) have been successfully
4 installed. If not, please consult the [Installation](@ref Installation)
5 section of this manual. If you are new to the [Julia](https://julialang.org)
6 programming language, the official manual has a useful guide to [getting
7 started with
8 Julia](https://docs.julialang.org/en/latest/manual/getting-started/).
9
10 In the following, two simple examples are presented using some of the core
11 functionality of Granular.jl. For more examples, see the scripts in the
12 [examples/](https://github.com/anders-dc/Granular.jl/tree/master/examples)
13 directory.
14
15 The relevant functions are all contained in the `Granular` module, which can be
16 imported with `import Granular`. *Note:* As per Julia conventions, functions
17 that contain an exclamation mark (!) modify the values of the arguments.
18
19 All of the functions called below are documented in the source code, and their
20 documentation can be found in the [Public API Index](@ref main-index) in the
21 online documentation, or simply from the Julia shell by typing `?<function
22 name>`. An example:
23
24 ```julia-repl
25 julia> ?Granular.createSimulation
26 createSimulation([id])
27
28 Create a simulation object to contain all relevant variables such as temporal
29 parameters, fluid grids, grains, and contacts. The parameter id is used to
30 uniquely identify the simulation when it is written to disk.
31
32 The function returns a Simulation object, which you can add grains to, e.g.
33 with addGrainCylindrical!.
34
35 Optional argument
36 ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
37
38 • id::String="unnamed": simulation identifying string.
39 ```
40
41 You can go through the examples below by typing in the lines starting with
42 `julia>` into the Julia interactive shell, which comes up when you start the
43 Julia app or run `julia` from the command line in a terminal. Do not include
44 the `julia> ` part, just the remaining text of that line.
45
46 Alternatively, you can create a Julia script with the file name extension
47 `.jl`. This file should contains all of the relevant commands in succession,
48 which is useful for quickly repeating runs. Julia scripts can be evaluated
49 form the command line using `julia <scriptname>.jl`.
50
51 ## Collision between two particles
52 For the first example (`examples/two-grains.jl`), we will create two grains and
53 make the first grain bump in to the second grain.
54
55 As the first step, we import all the Granular.jl functionality:
56
57 ```julia-repl
58 julia> import Granular
59 ```
60
61 ### Simulation setup
62 Next, we create our simulation object which holds all information related to
63 the simulated grains. In this documentation, we will use the name `sim` for
64 the simulation object:
65
66 ```julia-repl
67 julia> sim = Granular.createSimulation(id="two-grains")
68 Granular.Simulation("two-grains", 0, 0.0, -1.0, -1.0, -1.0, 0, 0.0,
69 Granular.GrainCylindrical[], Granular.Ocean(false, [0.0], [0.0], [0.0], [0.0],
70 [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], Array{Int64,1}[#undef], 1, 1,
71 1, 1), Granular.Atmosphere(false, [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
72 [0.0], [0.0], Array{Int64,1}[#undef], 1, 1, 1, 1, false), 16)
73 ```
74
75 After creating the simulation object `sim`, we will be presented with some
76 output about the default contents of the `sim` simulation object. This is of
77 minor importance as of now, and can safely be ignored.
78
79 In the above `createSimulation` call, the `id` argument is optional. It is used
80 to name simulation output files that are written to the disk.
81
82 ### Adding grains one by one
83 We have now created a simulation object `sim`, which will be used during all of
84 the following commands. Next, we can add grains to this object. The first
85 grain is cylinder shaped, placed at the x-y position (0,0) m. It has a radius
86 of 0.1 m, and has a thickness of 0.05 m. As this call modifies the `sim`
87 object, the function contains an exclamation mark (!). For further information
88 regarding this call, see the reference to [`addGrainCylindrical!`](@ref), found
89 in the [Public API documentation](@ref).
90
91 ```julia-repl
92 julia> Granular.addGrainCylindrical!(sim, [0.0, 0.0], 0.1, 0.05)
93 INFO: Added Grain 1
94 ```
95
96 Let's add another grain, placed at some distance from the first grain:
97
98 ```julia-repl
99 julia> Granular.addGrainCylindrical!(sim, [0.5, 0.0], 0.1, 0.05)
100 INFO: Added Grain 2
101 ```
102
103 We now want to prescribe a linear (not rotational or angular) velocity to the
104 first grain to make it bump into the second grain.
105
106 The simulation object `sim` contains an array of all grains that are added to
107 it. We can directly inspect the grains and their values from the simulation
108 object. Let's take a look at the default value of the linear velocity, called
109 `lin_vel`:
110
111 ```julia-repl
112 julia> sim.grains[1].lin_vel
113 2-element Array{Float64, 1}:
114 0.0
115 0.0
116 ```
117
118 The default value is a (0,0) vector, which means that it is not moving in
119 space. With a similar call, we can modify the properties of the first grain
120 directly and prescribe a velocity to it:
121
122 ```julia-repl
123 julia> sim.grains[1].lin_vel = [1.0, 0.0]
124 2-element Array{Float64, 1}:
125 1.0
126 0.0
127 ```
128
129 The first grain (index 1 in `sim.grains`) now has a positive velocity along `x`
130 with the value of 1.0 meter per second.
131
132 ### Setting temporal parameters for the simulation
133 Before we can start the simulation we need to set some more parameters vital to
134 the simulation, like what time step to use, how often to write output files to
135 the disk, and for how long to run the simulation. To set the computational
136 time step, we call the following:
137
138 ```julia-repl
139 julia> Granular.setTimeStep!(sim)
140 INFO: Time step length t=8.478741928617433e-5 s
141 ```
142
143 A suitable time step is automatically determined based on the size and elastic
144 properties of the grain.
145
146 The two grains are 0.3 meters apart, as the centers are placed 0.5 meter from
147 each other and each grain has a radius of 0.1 m. With a linear velocity of 1.0
148 m/s, the collision should occur after 0.3 seconds. For that reason it seems
149 reasonable to run the simulation for a total of 1.0 seconds. We choose to
150 produce an output file every 0.05 seconds, but this can be tweaked to taste.
151 We will later use the produced output files for visualization purposes.
152
153 ```julia-repl
154 julia> Granular.setOutputFileInterval!(sim, 0.05)
155
156 julia> Granular.setTotalTime!(sim, 1.0)
157 ```
158
159 ### Running the simulation
160 We are now ready to run the simulation. For illustrative purposes, let us
161 compare the kinetic energy in the granular system before and after the
162 collision. For now, we compute the initial value using the following call:
163
164 ```julia-repl
165 julia> Granular.totalGrainKineticTranslationalEnergy(sim)
166 0.7335618846132168
167 ```
168
169 The above value is the total translational (not angular) kinetic energy in
170 Joules before the simulation is started.
171
172 We have two choices for running the simulation; we can either run the entire
173 simulation length with a single call, which steps time until the total time is
174 reached and generates output files along the way. Alternatively, we can run
175 the simulation for a single time step a time, and inspect the progress or do
176 other modifications along the way.
177
178 Here, we will run the entire simulation in one go, and afterwards visualize the
179 grains from their output files using ParaView.
180
181 ```julia-repl
182 julia> Granular.run!(sim)
183
184 INFO: Output file: ./two-grains/two-grains.grains.1.vtu
185 INFO: wrote status to ./two-grains/two-grains.status.txt
186 t = 0.04239370964308682/1.0 s
187 INFO: Output file: ./two-grains/two-grains.grains.2.vtu
188 INFO: wrote status to ./two-grains/two-grains.status.txt
189
190 ...
191
192 INFO: Output file: ./two-grains/two-grains.grains.20.vtu
193 INFO: wrote status to ./two-grains/two-grains.status.txt
194 t = 0.9920128056483803/1.0 s
195 INFO: ./two-grains/two-grains.py written, execute with 'pvpython /Users/ad/code/Granular-ext/two-grains/two-grains.py'
196 INFO: wrote status to ./two-grains/two-grains.status.txt
197 t = 1.0000676104805686/1.0 s
198 ```
199
200 The code informs us of the simulation progress along the way. It also and
201 notifies us as output files are generated. This output can be disabled by
202 passing `verbose=false` to the `run!()` command. Finally, the code tells us
203 that it has generated a ParaView python file for visualization, called
204 `two-grains.py`, located in the `two-grains/` directory.
205
206 We are interested in getting an immediate idea of how the collision went,
207 before going further. We can print the new velocities with the following
208 commands:
209
210 ```julia-repl
211 julia> sim.grains[1].lin_vel
212 2-element Array{Float64, 1}:
213 7.58343e-5
214 0.0
215
216 julia> sim.grains[2].lin_vel
217 2-element Array{Float64, 1}:
218 0.999924
219 0.0
220 ```
221
222 The first grain has transferred effectively all of its kinetic energy to the
223 second grain during the cause of the simulation. The total kinetic energy now
224 is the following:
225
226 ```julia-repl
227 julia> Granular.totalGrainKineticTranslationalEnergy(sim)
228 0.7334506347624973
229 ```
230
231 The before and after values for total kinetic energy are reasonably close (to
232 less than 0.1 percent), which is what can be expected given the computation
233 accuracy of the algorithm.
234
235 ### Visualizing the output
236 To visualize the output we open [ParaView](https://www.paraview.org). The
237 output files of the simulation are written using the VTK (visualization
238 toolkit) format, which is natively supported by ParaView.
239
240 While the `.vtu` files produced during the simulation can be opened with
241 ParaView and visualized manually using *Glyph* filters, the simplest and
242 fastest way to visualize the data is to use the Python script generated for the
243 simulation by Granular.jl.
244
245 Open ParaView and open the *Python Shell*, found under the menu *Tools > Python
246 Shell*. In the pop-up dialog we select *Run Script*, which opens yet another
247 dialog prompting us to locate the visualization script
248 (`two-grains/two-grains.py` in our example). We locate this file, which is
249 placed under the directory from where we launched the `julia` session with the
250 commands above. If you are not sure what the current working directory for the
251 Julia session is, it can be displayed with the command `pwd()` in the Julia
252 shell.
253
254 After selecting the `two-grains/two-grains.py` script, we can close the *Python
255 Shell* window to inspect our simulation visually. Press the *Play* symbol in
256 the top toolbar, and see what happens.
257
258 By default, the grains are colored according to their size. Alternatively, you
259 can color the grains using different parameters, such as linear velocity,
260 number of contacts, etc. These parameters can be selected by changing the
261 chosen field under the *Glyph1* object in the *Pipeline Browser* on the left,
262 and by selecting a parameter for *Coloring*. Press the *Apply* button at the
263 top of the panel on the left to see the changes in effect.
264
265 **Tip:** If you have the command `pvpython` (ParaView Python) is available from
266 the command line, you can visualize the simulation directly from the Julia
267 shell without entering ParaView by the command `Granular.render(sim)`. The
268 program `pvpython` is included in the ParaView download, and is in the macOS
269 application bundle located in
270 `/Applications/Paraview-5.4.0.app/Contents/bin/pvpython`. Furthermore, if you
271 have the `convert` command from ImageMagick installed (`brew install
272 imagemagick` on macOS), the output images will be merged into an animated GIF.
273
274 ### Exercises
275 To gain more familiarity with the simulation procedure, I suggest experimenting
276 with the following exercises. *Tip:* You can rerun an experiment after
277 changing one or more parameters by calling `Granular.resetTime!(sim);
278 Granular.run!(sim)`. Changes in grain size, grain mass, or elastic properties
279 require a recomputation of the numerical time step
280 (`Granular.setTimeStep!(sim)`) before calling `Granular.run!(sim)`. The output
281 files will be overwritten, and changes will be immediately available in
282 ParaView.
283
284 - What effect does the grain size have on the time step?
285 - Try to make an oblique collision by placing one of the grains at a different
286 `y` position.
287 - What happens if the second grain is set to be fixed in space
288 (`sim.grains[2].fixed = true`)?
289 - How is the relationship between total kinetic energy before and after
290 affected by the choice of time step length? Try setting different time
291 step values, e.g. with `sim.time_step = 0.1234` and rerun the simulation.
292
293 ## Sedimentation of grains
294 Grains are known to settle under gravitation in water according to *Stoke's
295 law*, where resistive drag acts opposite of gravity and with a magnitude
296 according to the square root of velocity difference between water and grain.
297 Granular.jl offers simple fluid grids with prescribed velocity fields, and the
298 grains are met with drag in this grid.
299
300 In this example (`examples/sedimentation.jl`) we will initialize grains with a
301 range of grain sizes in a loose configuration, add gravity and a surrounding
302 fluid grid, and let the grains settle towards the bottom.
303
304 As in the previous example, we start by creating a fluid grid:
305
306 ```julia-repl
307 julia> import Granular
308 julia> sim = Granular.createSimulation(id="sedimentation")
309 ```
310
311 ### Creating a pseudo-random grain packing
312 Instead of manually adding grains one by one, we can use the
313 `regularPacking!()` function to add a regular grid of random-sized grains to
314 the simulation. Below, we specify that we want the grid of grains to be 7
315 grains wide along x, and 25 grains tall along y. We also specify the grain
316 radii to fall between 0.02 and 0.2 m. The sizes will be drawn from a power-law
317 distribution by default.
318
319 ```julia-repl
320 julia> Granular.regularPacking!(sim, [7, 25], 0.02, 0.2)
321 ```
322
323 Since we haven't explicitly set the grain sizes for this example, we can
324 inspect the values by plotting a histogram of sizes (only works if the `PyPlot`
325 package is installed with `Pkg.add("PyPlot")`):
326
327 ```julia-repl
328 julia> Granular.plotGrainSizeDistribution(sim)
329 INFO: sedimentation-grain-size-distribution.png
330 ```
331
332 The output informs us that we have the plot saved as an image with the file
333 name `sedimentation-grain-size-distribution.png`.
334
335 ### Creating a fluid grid
336 We can now create a fluid (ocean) grid spanning the extent of the grains
337 created above:
338
339 ```julia-repl
340 julia> Granular.fitGridToGrains!(sim, sim.ocean)
341 INFO: Created regular Granular.Ocean grid from [0.06382302477946442,
342 0.03387419706945263] to [3.0386621000253293, 10.87955941983313] with a cell
343 size of 0.3862075959573571 ([7, 28]).
344 ```
345
346 The code informs us of the number of grid cells in each dimension (7 by 28
347 cells), and the edge positions (x = 0.0638 to 3.04 m, y = 0.0339 to 10.9 m).
348
349 We want the boundaries of the above grid to be impermeable for the grains, so
350 they stack up at the bottom. Granular.jl acknowledges the boundary types with
351 a confirmation message:
352
353 ```julia-repl
354 julia> Granular.setGridBoundaryConditions!(sim.ocean, "impermeable")
355 West (-x): impermeable (3)
356 East (+x): impermeable (3)
357 South (-y): impermeable (3)
358 North (+y): impermeable (3)
359 ```
360
361 ### Adding gravitational acceleration
362 If we start the simulation now nothing would happen as gravity is disabled by
363 default. We can enable gravitational acceleration as a constant body force for
364 each grain (`Force = mass * acceleration`):
365
366 ```julia-repl
367 julia> g = [0.0, -9.8];
368 julia> for grain in sim.grains
369 Granular.addBodyForce!(grain, grain.mass*g)
370 end
371 ```
372
373 ### Setting temporal parameters
374 As before, we ask the code to select a suitable computational time step based
375 on grain sizes and properties:
376
377 ```julia-repl
378 julia> Granular.setTimeStep!(sim)
379 INFO: Time step length t=1.6995699879716792e-5 s
380 ```
381
382 We also, as before, set the total simulation time as well as the output file
383 interval:
384
385 ```julia-repl
386 julia> Granular.setTotalTime!(sim, 10.0)
387 julia> Granular.setOutputFileInterval!(sim, 0.2)
388 ```
389
390 ### Running the simulation
391 We are now ready to run the simulation:
392
393 ```julia-repl
394 julia> Granular.run!(sim)
395 INFO: Output file: ./sedimentation/sedimentation.grains.1.vtu
396 INFO: Output file: ./sedimentation/sedimentation.ocean.1.vts
397 INFO: wrote status to ./sedimentation/sedimentation.status.txt
398 t = 0.19884968859273294/10.0 s
399 INFO: Output file: ./sedimentation/sedimentation.grains.2.vtu
400 INFO: Output file: ./sedimentation/sedimentation.ocean.2.vts
401 INFO: wrote status to ./sedimentation/sedimentation.status.txt
402 t = 0.3993989471735396/10.0 s
403
404 ...
405
406 INFO: Output file: ./sedimentation/sedimentation.grains.50.vtu
407 INFO: Output file: ./sedimentation/sedimentation.ocean.50.vts
408 INFO: wrote status to ./sedimentation/sedimentation.status.txt
409 t = 9.998435334626701/10.0 s
410 INFO: ./sedimentation/sedimentation.py written, execute with 'pvpython /Users/ad/code/Granular-ext/examples/sedimentation/sedimentation.py'
411 INFO: wrote status to ./sedimentation/sedimentation.status.txt
412 t = 10.00001593471549/10.0 s
413 ```
414
415 The output can be plotted in ParaView as described in the `two-grain` example
416 above, or, if `pvpython` is available from the command line, directly from
417 Julia with the following command:
418
419 ```julia-repl
420 julia> Granular.render(sim, trim=false)
421 ```
422
423 ### Exercises
424 - How are the granular contact pressures distributed in the final result? You
425 can visualize this by selecting "Contact Pressure [Pa]" in the *Coloring*
426 field inside ParaView.
427 - Try running the above example, but without fluid drag. Disable the drag by
428 including the call `Granlar.disableOceanDrag!(grain)` in the `for` loop
429 where gravitational acceleration is set for each grain.
430 - How does the range of grain sizes affect the result? Try making all grains
431 bigger or smaller.
432 - How is the model performance effected if the grain-size distribution is
433 wide or narrow?
434 - Create a landslide by turning the gravitational acceleration vector (set the
435 `y` component to a non-zero value, and setting the side boundaries to be
436 periodic with `Granular.setGridBoundaryConditions!(sim.ocean, "periodic",
437 "east west")`.