Granular.jl

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

commit be609c80406b4b964b420d58d939fb3dd6e62647
parent 61a6dfbb5653a440e9a76183d34f059e065baf95
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Tue, 16 Jan 2018 12:50:31 -0800

Merge branch 'master' of github.com:anders-dc/Granular.jl

Diffstat:
MMakefile | 12+++++++++---
Msrc/datatypes.jl | 2+-
Msrc/grain.jl | 172+++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------
Msrc/io.jl | 28+++++++++++++++++++++++++---
Mtest/collision-2floes-normal.jl | 73+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mtest/runtests.jl | 6+++---
Mtest/vtk.jl | 2+-
7 files changed, 227 insertions(+), 68 deletions(-)

diff --git a/Makefile b/Makefile @@ -1,15 +1,21 @@ default: test .PHONY: test -test: test-julia-0.6 test-julia-0.7 +test: test-julia-0.6 #test-julia-0.7 .PHONY: test-julia-0.6 test-julia-0.6: - julia --color=yes -e 'Pkg.test("Granular")' + @#julia --color=yes -e 'Pkg.test("Granular")' + julia --color=yes -e 'Pkg.test("Granular")' \ + && notify-send Granular.jl tests completed successfully on Julia 0.6 \ + || notify-send Granular.jl failed on Julia 0.6 .PHONY: test-julia-0.7 test-julia-0.7: - julia-0.7 --color=yes -e 'Pkg.test("Granular")' + @#julia-0.7 --color=yes -e 'Pkg.test("Granular")' + julia-0.7 --color=yes -e 'Pkg.test("Granular")' \ + && notify-send Granular.jl tests completed successfully on Julia 0.7 \ + || notify-send Granular.jl failed on Julia 0.7 .PHONY: docs docs: diff --git a/src/datatypes.jl b/src/datatypes.jl @@ -138,7 +138,7 @@ mutable struct GrainArrays youngs_modulus::Vector{Float64} poissons_ratio::Vector{Float64} tensile_strength::Vector{Float64} - #tensile_heal_rate::Vector{Float64} + tensile_heal_rate::Vector{Float64} compressive_strength_prefactor::Vector{Float64} ocean_drag_coeff_vert::Vector{Float64} diff --git a/src/grain.jl b/src/grain.jl @@ -331,63 +331,119 @@ arrays in an `GrainArrays` structure. function convertGrainDataToArrays(simulation::Simulation) ifarr = GrainArrays( - Array{Float64}(length(simulation.grains)), - - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - - Array{Int}(length(simulation.grains)), - Array{Int}(length(simulation.grains)), - Array{Int}(length(simulation.grains)), - Array{Int}(length(simulation.grains)), - Array{Int}(length(simulation.grains)), - - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - Array{Float64}(length(simulation.grains)), - - Array{Float64}(length(simulation.grains)), - Array{Int}(length(simulation.grains)), - - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - zeros(Float64, 3, length(simulation.grains)), - - Array{Int}(length(simulation.grains)), - - ) + # Material properties + ## density + Array{Float64}(length(simulation.grains)), + + # Geometrical properties + ## thickness + Array{Float64}(length(simulation.grains)), + ## contact_radius + Array{Float64}(length(simulation.grains)), + ## areal_radius + Array{Float64}(length(simulation.grains)), + ## circumreference + Array{Float64}(length(simulation.grains)), + ## horizontal_surface_area + Array{Float64}(length(simulation.grains)), + ## side_surface_area + Array{Float64}(length(simulation.grains)), + ## volume + Array{Float64}(length(simulation.grains)), + ## mass + Array{Float64}(length(simulation.grains)), + ## moment_of_inertia + Array{Float64}(length(simulation.grains)), + + # Linear kinematic degrees of freedom along horiz plane + ## lin_pos + zeros(Float64, 3, length(simulation.grains)), + ## lin_vel + zeros(Float64, 3, length(simulation.grains)), + ## lin_acc + zeros(Float64, 3, length(simulation.grains)), + ## force + zeros(Float64, 3, length(simulation.grains)), + ## external_body_force + zeros(Float64, 3, length(simulation.grains)), + ## lin_disp + zeros(Float64, 3, length(simulation.grains)), + + # Angular kinematic degrees of freedom for vert. rot. + ## ang_pos + zeros(Float64, 3, length(simulation.grains)), + ## ang_vel + zeros(Float64, 3, length(simulation.grains)), + ## ang_acc + zeros(Float64, 3, length(simulation.grains)), + ## torque + zeros(Float64, 3, length(simulation.grains)), + + # Kinematic constraint flags + ## fixed + Array{Int}(length(simulation.grains)), + ## allow_x_acc + Array{Int}(length(simulation.grains)), + ## allow_y_acc + Array{Int}(length(simulation.grains)), + ## rotating + Array{Int}(length(simulation.grains)), + ## enabled + Array{Int}(length(simulation.grains)), + + # Rheological parameters + ## contact_stiffness_normal + Array{Float64}(length(simulation.grains)), + ## contact_stiffness_tangential + Array{Float64}(length(simulation.grains)), + ## contact_viscosity_normal + Array{Float64}(length(simulation.grains)), + ## contact_viscosity_tangential + Array{Float64}(length(simulation.grains)), + ## contact_static_friction + Array{Float64}(length(simulation.grains)), + ## contact_dynamic_friction + Array{Float64}(length(simulation.grains)), + + ## youngs_modulus + Array{Float64}(length(simulation.grains)), + ## poissons_ratio + Array{Float64}(length(simulation.grains)), + ## tensile_strength + Array{Float64}(length(simulation.grains)), + ## tensile_heal_rate + Array{Float64}(length(simulation.grains)), + ## compressive_strength_prefactor + Array{Float64}(length(simulation.grains)), + + # Ocean/atmosphere interaction parameters + ## ocean_drag_coeff_vert + Array{Float64}(length(simulation.grains)), + ## ocean_drag_coeff_horiz + Array{Float64}(length(simulation.grains)), + ## atmosphere_drag_coeff_vert + Array{Float64}(length(simulation.grains)), + ## atmosphere_drag_coeff_horiz + Array{Float64}(length(simulation.grains)), + + # Interaction + ## pressure + Array{Float64}(length(simulation.grains)), + ## n_contacts + Array{Int}(length(simulation.grains)), + + ## granular_stress + zeros(Float64, 3, length(simulation.grains)), + ## ocean_stress + zeros(Float64, 3, length(simulation.grains)), + ## atmosphere_stress + zeros(Float64, 3, length(simulation.grains)), + + # Visualization parameters + ## color + Array{Int}(length(simulation.grains)), + + ) # fill arrays for i=1:length(simulation.grains) @@ -439,6 +495,7 @@ function convertGrainDataToArrays(simulation::Simulation) ifarr.youngs_modulus[i] = simulation.grains[i].youngs_modulus ifarr.poissons_ratio[i] = simulation.grains[i].poissons_ratio ifarr.tensile_strength[i] = simulation.grains[i].tensile_strength + ifarr.tensile_heal_rate[i] = simulation.grains[i].tensile_heal_rate ifarr.compressive_strength_prefactor[i] = simulation.grains[i].compressive_strength_prefactor @@ -510,6 +567,7 @@ function deleteGrainArrays!(ifarr::GrainArrays) ifarr.youngs_modulus = f1 ifarr.poissons_ratio = f1 ifarr.tensile_strength = f1 + ifarr.tensile_heal_rate = f1 ifarr.compressive_strength_prefactor = f1 ifarr.ocean_drag_coeff_vert = f1 diff --git a/src/io.jl b/src/io.jl @@ -376,6 +376,10 @@ function writeGrainVTK(simulation::Simulation, WriteVTK.vtk_point_data(vtkfile, ifarr.torque, "Sum of torques [N*m]") WriteVTK.vtk_point_data(vtkfile, ifarr.fixed, "Fixed in space [-]") + WriteVTK.vtk_point_data(vtkfile, ifarr.allow_x_acc, + "Fixed but allow (x) acceleration [-]") + WriteVTK.vtk_point_data(vtkfile, ifarr.allow_y_acc, + "Fixed but allow (y) acceleration [-]") WriteVTK.vtk_point_data(vtkfile, ifarr.rotating, "Free to rotate [-]") WriteVTK.vtk_point_data(vtkfile, ifarr.enabled, "Enabled [-]") @@ -398,6 +402,8 @@ function writeGrainVTK(simulation::Simulation, "Poisson's ratio [-]") WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_strength, "Tensile strength [Pa]") + WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_heal_rate, + "Tensile healing rate [1/s]") WriteVTK.vtk_point_data(vtkfile, ifarr.compressive_strength_prefactor, "Compressive strength prefactor [m^0.5 Pa]") @@ -412,7 +418,6 @@ function writeGrainVTK(simulation::Simulation, WriteVTK.vtk_point_data(vtkfile, ifarr.pressure, "Contact pressure [Pa]") - WriteVTK.vtk_point_data(vtkfile, ifarr.n_contacts, "Number of contacts [-]") @@ -447,8 +452,8 @@ Saves grain interactions to `.vtp` files for visualization with VTK, for example in Paraview. Convert Cell Data to Point Data and use with Tube filter. """ function writeGrainInteractionVTK(simulation::Simulation, - filename::String; - verbose::Bool=false) + filename::String; + verbose::Bool=false) i1 = Int64[] i2 = Int64[] @@ -460,6 +465,7 @@ function writeGrainInteractionVTK(simulation::Simulation, tensile_stress = Float64[] shear_displacement = Vector{Float64}[] contact_age = Float64[] + for i=1:length(simulation.grains) for ic=1:simulation.Nc_max if simulation.grains[i].contacts[ic] > 0 @@ -476,6 +482,19 @@ function writeGrainInteractionVTK(simulation::Simulation, r_i = simulation.grains[i].contact_radius r_j = simulation.grains[j].contact_radius + + # skip visualization of contacts across periodic BCs + if dist > 5.0*(r_i + r_j) && + (simulation.ocean.bc_west == 2 || + simulation.ocean.bc_east == 2 || + simulation.ocean.bc_north == 2 || + simulation.ocean.bc_south == 2 || + simulation.atmosphere.bc_west == 2 || + simulation.atmosphere.bc_east == 2 || + simulation.atmosphere.bc_north == 2 || + simulation.atmosphere.bc_south == 2) + continue + end δ_n = dist - (r_i + r_j) R_ij = harmonicMean(r_i, r_j) @@ -809,6 +828,8 @@ imagegrains.PointArrayStatus = [ 'Angular acceleration [rad s^-2]', 'Sum of torques [N*m]', 'Fixed in space [-]', +'Fixed but allow (x) acceleration [-]', +'Fixed but allow (y) acceleration [-]', 'Free to rotate [-]', 'Enabled [-]', 'Contact stiffness (normal) [N m^-1]', @@ -820,6 +841,7 @@ imagegrains.PointArrayStatus = [ "Young's modulus [Pa]", "Poisson's ratio [-]", 'Tensile strength [Pa]' +'Tensile heal rate [1/s]' 'Compressive strength prefactor [m^0.5 Pa]', 'Ocean drag coefficient (vertical) [-]', 'Ocean drag coefficient (horizontal) [-]', diff --git a/test/collision-2floes-normal.jl b/test/collision-2floes-normal.jl @@ -273,3 +273,76 @@ E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) @test E_kin_rot_init ≈ E_kin_rot_final @test sim.grains[2].lin_pos[1] > grain2_pos_init[1] +r = 10. +for angle in linspace(0, 2π, 4) + info("## Contact angle: $angle") + + info("Testing behavior with two fixed grains and allow_*_acc") + sim = Granular.createSimulation(id="test") + Granular.addGrainCylindrical!(sim, [0., 0.], r, 1., verbose=verbose) + Granular.addGrainCylindrical!(sim, [2.0*r*cos(angle), 2.0*r*sin(angle)], + r, 1., verbose=verbose) + sim.grains[1].fixed = true + sim.grains[2].fixed = true + + E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) + E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) + grain1_pos_init = sim.grains[1].lin_pos + grain2_pos_init = sim.grains[2].lin_pos + + Granular.setTotalTime!(sim, 10.0) + Granular.setTimeStep!(sim, epsilon=0.07) + sim_init = deepcopy(sim) + + info("Both fixed, no allow_*_acc, no cohesion (TY2)") + sim = deepcopy(sim_init) + #sim.grains[2].allow_y_acc = true # should not influence result + tol = 0.02 + Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) + E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) + E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) + @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol + @test E_kin_rot_init ≈ E_kin_rot_final + @test sim.grains[1].lin_pos ≈ grain1_pos_init + @test sim.grains[2].lin_pos ≈ grain2_pos_init + + info("Both fixed, no allow_*_acc, no cohesion (TY3)") + sim = deepcopy(sim_init) + #sim.grains[2].allow_y_acc = true # should not influence result + tol = 0.02 + Granular.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verbose) + E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) + E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) + @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol + @test E_kin_rot_init ≈ E_kin_rot_final + @test sim.grains[1].lin_pos ≈ grain1_pos_init + @test sim.grains[2].lin_pos ≈ grain2_pos_init + + info("Both fixed, no allow_*_acc, cohesion (TY2)") + sim = deepcopy(sim_init) + #sim.grains[2].allow_y_acc = true # should not influence result + sim.grains[1].tensile_strength = 200e3 + sim.grains[2].tensile_strength = 200e3 + tol = 0.2 + Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) + E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) + E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) + @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol + @test E_kin_rot_init ≈ E_kin_rot_final + @test sim.grains[1].lin_pos ≈ grain1_pos_init + @test sim.grains[2].lin_pos ≈ grain2_pos_init + + info("Both fixed, no allow_*_acc, cohesion (TY3)") + sim = deepcopy(sim_init) + #sim.grains[2].allow_y_acc = true # should not influence result + sim.grains[1].tensile_strength = 200e3 + sim.grains[2].tensile_strength = 200e3 + tol = 0.02 + Granular.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verbose) + E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) + E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) + @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol + @test E_kin_rot_init ≈ E_kin_rot_final + @test sim.grains[1].lin_pos ≈ grain1_pos_init + @test sim.grains[2].lin_pos ≈ grain2_pos_init +end diff --git a/test/runtests.jl b/test/runtests.jl @@ -1,6 +1,9 @@ using Compat.Test import Granular +include("collision-2floes-normal.jl") +include("collision-5floes-normal.jl") +include("collision-2floes-oblique.jl") include("grid.jl") include("contact-search-and-geometry.jl") include("grid-boundaries.jl") @@ -11,9 +14,6 @@ include("grain.jl") include("packing.jl") include("util.jl") include("temporal.jl") -include("collision-2floes-normal.jl") -include("collision-5floes-normal.jl") -include("collision-2floes-oblique.jl") include("cohesion.jl") include("netcdf.jl") include("vtk.jl") diff --git a/test/vtk.jl b/test/vtk.jl @@ -29,7 +29,7 @@ end grainpath = "test/test.grains.1.vtu" grainchecksum = -"c1e85bf2c5379050ba6dd674eecedaf6857f309015e94544d1dc27e8b46607c8 " * +"f5a5106bffaaa2fd556ab0cac602127aecdb9cd9fc74141b75d1a7e7f39e21d6 " * grainpath * "\n" graininteractionpath = "test/test.grain-interaction.1.vtp"