Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit aa49a17a7ce969f67b58e7f4196f65a48a5f0e5c
parent 5b777d9bdb78d24eab9387b0ab15e99ba074737f
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Thu, 20 Apr 2017 15:51:52 -0400

add docstrings to all functions

Diffstat:
Msrc/SeaIce.jl | 25+++++--------------------
Dsrc/arrays.jl | 1-
Msrc/contact_search.jl | 7+++++++
Msrc/datatypes.jl | 1-
Msrc/interaction.jl | 4++++
Msrc/io.jl | 5+++++
Msrc/simulation.jl | 67+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------
Msrc/temporal.jl | 27+++++++++++++++++++++++++--
Msrc/temporal_integration.jl | 18+++++++++++++++++-
Mtest/contact-search-and-geometry.jl | 4++--
10 files changed, 126 insertions(+), 33 deletions(-)

diff --git a/src/SeaIce.jl b/src/SeaIce.jl @@ -1,30 +1,15 @@ #!/usr/bin/env julia -############################## -# Sea-ice dynamics simulator # -############################## - -#= -Add SeaIce to your julia path, e.g. by: - push!(LOAD_PATH, "/home/user/src/seaice/") - -If this statement is added to `~/.juliarc.jl`, it will become persistent between -julia sessions. Note that the `~` symbol for the home folder does not seem to -work (julia v. 0.4.1) in the `.juliarc.jl` file. - -Import package contents with: - import SeaIce - -If this file is changed, reimport it with: - reload("SeaIce") -=# - +""" +# SeaIce.jl +Offline sea-ice dynamics simulator by [Anders Damsgaard](mailto:andersd@riseup.net), +[www.adamsgaard.dk](https://adamsgaard.dk) +""" module SeaIce include("datatypes.jl") include("icefloe.jl") include("simulation.jl") -include("arrays.jl") include("grid.jl") include("packing.jl") include("contact_search.jl") diff --git a/src/arrays.jl b/src/arrays.jl @@ -1 +0,0 @@ -## Global arrays for general simulation data diff --git a/src/contact_search.jl b/src/contact_search.jl @@ -1,4 +1,11 @@ ## Contact mapping +""" +Top-level function to perform an inter-ice floe contact search, based on ice +floe linear positions and contact radii. + +The simplest contact search algorithm (`method="all to all"`) is the most +computationally expensive (O(n^2)). +""" function findContacts!(simulation::Simulation, method::String = "all to all") diff --git a/src/datatypes.jl b/src/datatypes.jl @@ -3,7 +3,6 @@ const float = Float64 const vector = Array{Float64, 1} ## Particle composite types - type IceFloeCylindrical # Material properties diff --git a/src/interaction.jl b/src/interaction.jl @@ -45,6 +45,10 @@ function interactIceFloes!(simulation::Simulation, end +""" +Resolves linear-elastic interaction between two ice floes in the contact-normal +direction. +""" function interactNormalLinearElastic(simulation::Simulation, i::Integer, j::Integer, overlap_vector::vector) diff --git a/src/io.jl b/src/io.jl @@ -2,6 +2,11 @@ import WriteVTK # Install with Pkg.add("WriteVTK") ## IO functions +""" +Write a VTK file to disk containing all ice floes in the `simulation` in an +unstructured mesh (file type `.vtu`). These files can be read by ParaView and +can be visualized by applying a *Glyph* filter. +""" function writeVTK(simulation::Simulation; folder::String=".", verbose::Bool=false) diff --git a/src/simulation.jl b/src/simulation.jl @@ -1,5 +1,24 @@ ## General simulation functions +""" + createSimulation([id::String="unnamed", + time_iteration::Int=0, + time::Float64=0.0, + time_total::Float64=-1., + time_step::Float64=-1., + file_time_step::Float64=-1., + file_number::Int=0, + ice_floes=Array{IceFloeCylindrical, 1}[], + contact_pairs=Array{Int64, 1}[], + overlaps=Array{Array{Float64, 1}, 1}[], + wall_contacts=Array{Int64, 1}[]]) + +Create a simulation object containing all relevant variables such as temporal +parameters, and lists of ice floes and contacts. + +The parameter `id` is used to uniquely identify the simulation when it is +written to disk. +""" function createSimulation(;id::String="unnamed", time_iteration::Int=0, time::Float64=0.0, @@ -25,12 +44,43 @@ function createSimulation(;id::String="unnamed", wall_contacts) end +""" + run!(simulation[, + verbose::Bool = true, + status_interval = 100., + show_file_output = true, + single_step=false, + temporal_integration_method="Three-term Taylor"]) + +Run the `simulation` through time until `simulation.time` equals or exceeds +`simulatim.time_total`. This function requires that all ice floes are added to +the simulation and that the length of the computational time step is adjusted +accordingly. + +The function will search for contacts, determine the force balance on each ice +floe, and integrate all kinematic degrees of freedom accordingly. The temporal +integration is explicit and of length `simulation.time_step`. This function +will write VTK files to disk in the intervals `simulation.file_time_step` by the +function `writeVTK`. If this value is negative, no output files will be written +to disk. + +# Arguments +* `simulation::Simulation`: the simulation to run (object is modified) +* `verbose::Bool=true`: show verbose information during the time loop +* `status_interval::Bool=true`: show verbose information during the time loop +* `show_file_output::Bool=true`: report to stdout when output file is written +* `single_step::Bool=false`: run simulation for a single time step only. If + this causes `simulation.time` to exceed `simulation.time_total`, the latter + is increased accordingly. +* `temporal_integration_method::String="Three-term Taylor"`: type of integration + method to use. See `updateIceFloeKinematics` for details. +""" function run!(simulation::Simulation; - verbose::Bool = true, - status_interval = 100., - show_file_output = true, - single_step=false, - temporal_integration_method="Three-term Taylor") + verbose::Bool=true, + status_interval::Int=100., + show_file_output::Bool=true, + single_step::Bool=false, + temporal_integration_method::String="Three-term Taylor") if single_step && simulation.time >= simulation.time_total simulation.time_total += simulation.time_step @@ -63,7 +113,7 @@ function run!(simulation::Simulation; # Update time variables simulation.time_iteration += 1 - simulation.time += simulation.time_step + incrementCurrentTime(simulation, simulation.time_step) time_since_output_file = time_since_output_file + simulation.time_step if single_step @@ -79,6 +129,8 @@ function run!(simulation::Simulation; end end +"Add an `icefloe` to the `simulation` object. If `verbose` is true, a short +confirmation message will be printed to stdout`." function addIceFloe!(simulation::Simulation, icefloe::IceFloeCylindrical, verbose::Bool = False) @@ -90,6 +142,7 @@ function addIceFloe!(simulation::Simulation, end end +"Remove ice floe with index `i` from the `simulation` object." function removeIceFloe!(simulation::Simulation, i::Integer) if i < 1 error("Index must be greater than 0 (i = $i)") @@ -98,6 +151,7 @@ function removeIceFloe!(simulation::Simulation, i::Integer) delete!(simulation.ice_floes, i) end +"Sets the `force` and `torque` values of all ice floes to zero." function zeroForcesAndTorques!(simulation::Simulation) for icefloe in simulation.ice_floes icefloe.force = zeros(2) @@ -105,6 +159,7 @@ function zeroForcesAndTorques!(simulation::Simulation) end end +"Prints the current simulation time and total time to standard out" function reportSimulationTimeToStdout(simulation::Simulation) print("\r t = ", simulation.time, '/', simulation.time_total, " s ") diff --git a/src/temporal.jl b/src/temporal.jl @@ -1,4 +1,9 @@ -## Set temporal parameters +""" + setTotalTime!(simulation::Simulation, t::float) + +Sets the total simulation time of the `simulation` object to `t`, with parameter +value checks. +""" function setTotalTime!(simulation::Simulation, t::float) if t <= 0.0 error("Total time should be a positive value (t = $t s)") @@ -6,6 +11,12 @@ function setTotalTime!(simulation::Simulation, t::float) simulation.time_total = t end +""" + setCurrentTime!(simulation::Simulation, t::float) + +Sets the current simulation time of the `simulation` object to `t`, with +parameter value checks. +""" function setCurrentTime!(simulation::Simulation, t::float) if t <= 0.0 error("Current time should be a positive value (t = $t s)") @@ -13,11 +24,17 @@ function setCurrentTime!(simulation::Simulation, t::float) simulation.time = t end +""" + incrementCurrentTime!(simulation::Simulation, t::float) + +Sets the current simulation time of the `simulation` object to `t`, with +parameter value checks. +""" function incrementCurrentTime!(simulation::Simulation, t::float) if t <= 0.0 error("Current time increment should be a positive value (t = $t s)") end - simulation.time = simulation.time + t + simulation.time += t end function setOutputFileInterval!(simulation::Simulation, t::float) @@ -27,10 +44,12 @@ function setOutputFileInterval!(simulation::Simulation, t::float) simulation.file_time_step = t end +"Disables the write of output files to disk during a simulation." function disableOutputFiles!(simulation::Simulation) simulation.file_time_step = 0.0 end +"Checks if simulation parameters are of reasonable values." function checkTimeParameters(simulation::Simulation) if simulation.time_total <= 0.0 || simulation.time_total <= simulation.time error("Total time should be positive and larger than current time.\n", @@ -42,6 +61,8 @@ function checkTimeParameters(simulation::Simulation) end end +"Finds the smallest mass of all ice floes in a simulation. Used to determine +the optimal time step length." function findSmallestIceFloeMass(simulation::Simulation) m_min = 1e20 i_min = -1 @@ -55,6 +76,8 @@ function findSmallestIceFloeMass(simulation::Simulation) return m_min, i_min end +"Finds the largest elastic stiffness of all ice floes in a simulation. Used to +determine the optimal time step length." function findLargestIceFloeStiffness(simulation::Simulation) k_n_max = 0. k_t_max = 0. diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl @@ -1,6 +1,19 @@ """ + updateIceFloeKinematics!(simulation::Simulation[, + method::String = "Three-term Taylor"]) + Update the ice floe kinematic parameters using a temporal integration scheme, the current force and torque balance, and gravitational acceleration. + +# Arguments +* `simulation::Simulation`: update the ice floe positions in this object + according to temporal integration of length `simulation.time_step`. +* `method::String = "Three-term Taylor"`: the integration method to use. + Available methods are "Two-term Taylor" and "Three-term Taylor". The + three-term Taylor expansion is slightly more computationally expensive than + the two-term Taylor expansion, but offers an order-of-magnitude increase in + precision of ice floe positions. The two-term expansion can obtain similar + precision if the time step is 1/10 the length. """ function updateIceFloeKinematics!(simulation::Simulation; method::String = "Three-term Taylor") @@ -18,6 +31,8 @@ function updateIceFloeKinematics!(simulation::Simulation; end end +"Use a two-term Taylor expansion for integrating the kinematic degrees of +freedom for an `icefloe`." function updateIceFloeKinematicsTwoTermTaylor(icefloe::IceFloeCylindrical, simulation::Simulation) icefloe.lin_acc = icefloe.force/icefloe.mass @@ -41,6 +56,8 @@ function updateIceFloeKinematicsTwoTermTaylor(icefloe::IceFloeCylindrical, icefloe.ang_vel += icefloe.ang_acc * simulation.time_step end +"Use a three-term Taylor expansion for integrating the kinematic degrees of +freedom for an `icefloe`." function updateIceFloeKinematicsThreeTermTaylor(icefloe::IceFloeCylindrical, simulation::Simulation) @@ -80,4 +97,3 @@ function updateIceFloeKinematicsThreeTermTaylor(icefloe::IceFloeCylindrical, icefloe.ang_vel += icefloe.ang_acc * simulation.time_step + 0.5 * d_ang_acc_dt * simulation.time_step^2. end - diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl @@ -14,8 +14,8 @@ SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false) position_ij = SeaIce.interIceFloePositionVector(sim, 1, 2) overlap_ij = SeaIce.findOverlap(sim, 1, 2, position_ij) -Base.Test.@test_approx_eq [18., 0.] position_ij -Base.Test.@test_approx_eq -2. overlap_ij +Base.Test.@test [18., 0.] \approx position_ij +Base.Test.@test -2. \approx overlap_ij info("Testing findContactsAllToAll(...)")