tsimulation.jl - 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
---
tsimulation.jl (13648B)
---
1 using Printf
2
3 ## General simulation functions
4
5 export createSimulation
6 """
7 createSimulation([id])
8
9 Create a simulation object to contain all relevant variables such as temporal
10 parameters, fluid grids, grains, and contacts. The parameter `id` is used to
11 uniquely identify the simulation when it is written to disk.
12
13 The function returns a `Simulation` object, which you can add grains to, e.g.
14 with [`addGrainCylindrical!`](@ref).
15
16 # Optional argument
17 * `id::String="unnamed"`: simulation identifying string.
18
19 """
20 function createSimulation(;id::String="unnamed")
21
22 # default values for simulation parameters (not included as arguments as
23 # they are very rarely chagned and make the docstring much more cluttered
24 # and intimidating
25 time_iteration::Int = 0
26 time::Float64 = 0.0
27 time_total::Float64 = -1.
28 time_step::Float64 = -1.
29 file_time_step::Float64 = -1.
30 file_number::Int = 0
31 file_time_since_output_file::Float64 = 0.
32 grains = Array{GrainCylindrical, 1}[]
33 ocean::Ocean = createEmptyOcean()
34 atmosphere::Atmosphere = createEmptyAtmosphere()
35 Nc_max::Int = 32
36 walls = Array{WallLinearFrictionless, 1}[]
37
38 return Simulation(id,
39 time_iteration,
40 time,
41 time_total,
42 time_step,
43 file_time_step,
44 file_number,
45 file_time_since_output_file,
46 grains,
47 ocean,
48 atmosphere,
49 Nc_max,
50 walls)
51 end
52 function createSimulation(id::String)
53 createSimulation(id=id)
54 end
55
56 export run!
57 """
58 run!(simulation[,
59 verbose::Bool = true,
60 status_interval = 100.,
61 show_file_output = true,
62 single_step = false,
63 temporal_integration_method = "Three-term Taylor"],
64 write_jld2 = false)
65
66 Run the `simulation` through time until `simulation.time` equals or exceeds
67 `simulatim.time_total`. This function requires that all grains are added to
68 the simulation and that the length of the computational time step is adjusted
69 accordingly.
70
71 The function will search for contacts, determine the force balance on each ice
72 floe, and integrate all kinematic degrees of freedom accordingly. The temporal
73 integration is explicit and of length `simulation.time_step`. This function
74 will write VTK files to disk in the intervals `simulation.file_time_step` by the
75 function `writeVTK`. If this value is negative, no output files will be written
76 to disk.
77
78 # Arguments
79 * `simulation::Simulation`: the simulation to run (object is modified)
80 * `verbose::Bool=true`: show verbose information during the time loop
81 * `status_interval::Bool=true`: show verbose information during the time loop
82 * `show_file_output::Bool=true`: report to stdout when output file is written
83 * `single_step::Bool=false`: run simulation for a single time step only. If
84 this causes `simulation.time` to exceed `simulation.time_total`, the latter
85 is increased accordingly.
86 * `temporal_integration_method::String="Three-term Taylor"`: type of integration
87 method to use. See `updateGrainKinematics` for details.
88 * `write_jld2::Bool=false`: write simulation state to disk as JLD2 files (see
89 `Granular.writeSimulation(...)` whenever saving VTK output.
90 """
91 function run!(simulation::Simulation;
92 verbose::Bool=true,
93 status_interval::Int=100,
94 show_file_output::Bool=true,
95 single_step::Bool=false,
96 temporal_integration_method::String="Three-term Taylor",
97 write_jld2::Bool=false)
98
99 if single_step && simulation.time >= simulation.time_total
100 simulation.time_total += simulation.time_step
101 end
102
103 checkTimeParameters(simulation, single_step=single_step)
104
105 # if both are enabled, check if the atmosphere grid spatial geometry is
106 # identical to the ocean grid
107 if simulation.time_iteration == 0 &&
108 typeof(simulation.atmosphere.input_file) != Bool &&
109 typeof(simulation.ocean.input_file) != Bool
110
111 if simulation.ocean.xq ≈ simulation.atmosphere.xq &&
112 simulation.ocean.yq ≈ simulation.atmosphere.yq
113 if verbose
114 @info "identical ocean and atmosphere grids, " *
115 "turning on grid optimizations"
116 end
117 simulation.atmosphere.collocated_with_ocean_grid = true
118 end
119 end
120
121
122 # number of time steps between output files
123 n_file_time_step = Int(ceil(simulation.file_time_step/simulation.time_step))
124
125 while simulation.time <= simulation.time_total
126
127 if simulation.file_time_step > 0.0 &&
128 simulation.time_iteration % n_file_time_step == 0
129
130 if show_file_output
131 println()
132 end
133 if write_jld2
134 writeSimulation(simulation, verbose=show_file_output)
135 end
136 writeVTK(simulation, verbose=show_file_output)
137 writeSimulationStatus(simulation, verbose=show_file_output)
138 simulation.file_time_since_output_file = 0.0
139 end
140
141 if verbose && simulation.time_iteration % status_interval == 0
142 reportSimulationTimeToStdout(simulation)
143 end
144
145 zeroForcesAndTorques!(simulation)
146
147 if typeof(simulation.atmosphere.input_file) != Bool &&
148 !simulation.atmosphere.collocated_with_ocean_grid
149 sortGrainsInGrid!(simulation, simulation.atmosphere)
150 end
151
152 if typeof(simulation.ocean.input_file) != Bool
153 sortGrainsInGrid!(simulation, simulation.ocean)
154 findContacts!(simulation, method="ocean grid")
155
156 if simulation.atmosphere.collocated_with_ocean_grid
157 copyGridSortingInfo!(simulation.ocean, simulation.atmosphere,
158 simulation.grains)
159 end
160
161 elseif typeof(simulation.atmosphere.input_file) != Bool
162 findContacts!(simulation, method="atmosphere grid")
163
164 else
165 findContacts!(simulation, method="all to all")
166 end
167
168 interact!(simulation)
169 interactWalls!(simulation)
170
171 if typeof(simulation.ocean.input_file) != Bool
172 addOceanDrag!(simulation)
173 end
174
175 if typeof(simulation.atmosphere.input_file) != Bool
176 addAtmosphereDrag!(simulation)
177 end
178
179 updateGrainKinematics!(simulation, method=temporal_integration_method)
180 updateWallKinematics!(simulation, method=temporal_integration_method)
181
182 # Update time variables
183 simulation.time_iteration += 1
184 incrementCurrentTime!(simulation, simulation.time_step)
185
186 if single_step
187 return nothing
188 end
189 end
190
191 if simulation.file_time_step > 0.0
192 if show_file_output
193 println()
194 end
195 writeParaviewPythonScript(simulation, verbose=show_file_output)
196 writeSimulationStatus(simulation, verbose=show_file_output)
197 end
198
199 if verbose
200 reportSimulationTimeToStdout(simulation)
201 println()
202 end
203 nothing
204 end
205
206 export addGrain!
207 """
208 addGrain!(simulation::Simulation,
209 grain::GrainCylindrical,
210 verbose::Bool = false)
211
212 Add an `grain` to the `simulation` object. If `verbose` is true, a short
213 confirmation message will be printed to stdout.
214 """
215 function addGrain!(simulation::Simulation,
216 grain::GrainCylindrical,
217 verbose::Bool = false)
218 push!(simulation.grains, grain)
219
220 if verbose
221 @info "Added grain $(length(simulation.grains))"
222 end
223 nothing
224 end
225
226 export addWall!
227 """
228 addWall!(simulation::Simulation,
229 wall::WallLinearFrictionless,
230 verbose::Bool = false)
231
232 Add an `wall` to the `simulation` object. If `verbose` is true, a short
233 confirmation message will be printed to stdout.
234 """
235 function addWall!(simulation::Simulation,
236 wall::WallLinearFrictionless,
237 verbose::Bool = false)
238 push!(simulation.walls, wall)
239
240 if verbose
241 @info "Added wall $(length(simulation.walls))"
242 end
243 nothing
244 end
245
246 export disableGrain!
247 "Disable grain with index `i` in the `simulation` object."
248 function disableGrain!(simulation::Simulation, i::Int)
249 if i < 1
250 error("Index must be greater than 0 (i = $i)")
251 end
252 simulation.grains[i].enabled = false
253 nothing
254 end
255
256 export zeroForcesAndTorques!
257 "Sets the `force` and `torque` values of all grains and dynamic walls to zero."
258 function zeroForcesAndTorques!(simulation::Simulation)
259 for grain in simulation.grains
260 grain.force .= grain.external_body_force
261 grain.torque = zeros(3)
262 grain.pressure = 0.
263 end
264 for wall in simulation.walls
265 wall.force = 0.
266 end
267 nothing
268 end
269
270 export reportSimulationTimeToStdout
271 "Prints the current simulation time and total time to standard out"
272 function reportSimulationTimeToStdout(simulation::Simulation)
273 print("\r t = ", simulation.time, '/', simulation.time_total,
274 " s ")
275 nothing
276 end
277
278 export compareSimulations
279 """
280 compareSimulations(sim1::Simulation, sim2::Simulation)
281
282 Compare values of two `Simulation` objects using the `Base.Test` framework.
283 """
284 function compareSimulations(sim1::Simulation, sim2::Simulation)
285
286 Test.@test sim1.id == sim2.id
287
288 Test.@test sim1.time_iteration == sim2.time_iteration
289 Test.@test sim1.time ≈ sim2.time
290 Test.@test sim1.time_total ≈ sim2.time_total
291 Test.@test sim1.time_step ≈ sim2.time_step
292 Test.@test sim1.file_time_step ≈ sim2.file_time_step
293 Test.@test sim1.file_number == sim2.file_number
294 Test.@test sim1.file_time_since_output_file ≈ sim2.file_time_since_output_file
295
296 for i=1:length(sim1.grains)
297 compareGrains(sim1.grains[i], sim2.grains[i])
298 end
299 compareOceans(sim1.ocean, sim2.ocean)
300 compareAtmospheres(sim1.atmosphere, sim2.atmosphere)
301
302 Test.@test sim1.Nc_max == sim2.Nc_max
303 nothing
304 end
305
306 export printMemoryUsage
307 """
308 printMemoryUsage(sim::Simulation)
309
310 Shows the memory footprint of the simulation object.
311 """
312 function printMemoryUsage(sim::Simulation)
313
314 reportMemory(sim, "sim")
315 println(" where:")
316
317 reportMemory(sim.grains, " sim.grains",
318 "(N=$(length(sim.grains)))")
319
320 reportMemory(sim.walls, " sim.walls",
321 "(N=$(length(sim.walls)))")
322
323 reportMemory(sim.ocean, " sim.ocean",
324 "($(size(sim.ocean.xh, 1))x" *
325 "$(size(sim.ocean.xh, 2))x" *
326 "$(size(sim.ocean.xh, 3)))")
327
328 reportMemory(sim.atmosphere, " sim.atmosphere",
329 "($(size(sim.atmosphere.xh, 1))x" *
330 "$(size(sim.atmosphere.xh, 2))x" *
331 "$(size(sim.atmosphere.xh, 3)))")
332 nothing
333 end
334
335 function reportMemory(variable, head::String, tail::String="")
336 bytes = Base.summarysize(variable)
337 if bytes < 10_000
338 size_str = @sprintf "%5d bytes" bytes
339 elseif bytes < 10_000 * 1024
340 size_str = @sprintf "%5d kb" bytes ÷ 1024
341 elseif bytes < 10_000 * 1024 * 1024
342 size_str = @sprintf "%5d Mb" bytes ÷ 1024 ÷ 1024
343 else
344 size_str = @sprintf "%5d Gb" bytes ÷ 1024 ÷ 1024 ÷ 1024
345 end
346 @printf("%-20s %s %s\n", head, size_str, tail)
347 nothing
348 end
349
350 export setMaximumNumberOfContactsPerGrain!
351 """
352 setMaximumNumberOfContactsPerGrain!(simulation, number_of_contacts)
353
354 Change the maximum number of contacts per grain, which changes simulation.Nc_max
355 and reallocates memory for each grain. Larger values require more memory, but
356 allow simulation of wider grain-size distributions. The default value is a
357 maximum of 32 contacts per grain, which is sufficient for most practical
358 purposes.
359
360 # Arguments
361 * `simulation::Simulation`: the Simulation object to modify
362 * `number_of_contacts::Int`: the maximum number of contacts per grain to allow.
363 """
364 function setMaximumNumberOfContactsPerGrain!(sim::Simulation,
365 number_of_contacts::Int)
366
367 if number_of_contacts < 1
368 error("the parameter number_of_contacts must be a positive integer, " *
369 "but has the value '$number_of_contacts'")
370 end
371 if number_of_contacts == sim.Nc_max
372 error("number_of_contacts equals the current number of contacts " *
373 "sim.Nc_max = $(sim.Nc_max)")
374 end
375
376 Nc_max_orig = sim.Nc_max
377 sim.Nc_max = number_of_contacts
378 diff = sim.Nc_max - Nc_max_orig
379
380 for grain in sim.grains
381
382 if diff > 0
383 # push values to the end of contact arrays if Nc_max > Nc_max_orig
384 for i=1:diff
385 push!(grain.contacts, 0)
386 push!(grain.position_vector, zeros(Float64, 3))
387 push!(grain.contact_parallel_displacement, zeros(Float64, 3))
388 push!(grain.contact_rotation, zeros(Float64, 3))
389 push!(grain.contact_age, 0.0)
390 push!(grain.contact_area, 0.0)
391 push!(grain.compressive_failure, false)
392 end
393
394 else
395 # pop values from the end of contact arrays if Nc_max < Nc_max_orig
396 for i=1:abs(diff)
397 pop!(grain.contacts)
398 pop!(grain.position_vector)
399 pop!(grain.contact_parallel_displacement)
400 pop!(grain.contact_rotation)
401 pop!(grain.contact_age)
402 pop!(grain.contact_area)
403 pop!(grain.compressive_failure)
404 end
405 end
406 end
407 end