Granular.jl

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

commit 34985c227ef815ae8c6d65e0ceceaec7e69009d4
parent a25171d98558a1f674f705693935f7d75927ca33
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Thu,  9 Nov 2017 10:49:37 -0600

fix several syntax changes for julia 0.7

Diffstat:
Msrc/atmosphere.jl | 10+++++-----
Msrc/grid.jl | 13+++++++------
Msrc/interaction.jl | 4++--
Msrc/io.jl | 3++-
Msrc/ocean.jl | 10+++++-----
Msrc/packing.jl | 10+++++-----
Msrc/temporal.jl | 4++--
Msrc/temporal_integration.jl | 4++--
Msrc/util.jl | 7++++---
Mtest/netcdf.jl | 12++++++------
10 files changed, 40 insertions(+), 37 deletions(-)

diff --git a/src/atmosphere.jl b/src/atmosphere.jl @@ -220,11 +220,11 @@ function applyAtmosphereVorticityToGrain!(grain::GrainCylindrical, rho_a = 1.2754 # atmosphere density grain.torque += - pi*grain.areal_radius^4.*rho_a* - (grain.areal_radius/5.*grain.atmosphere_drag_coeff_horiz + - .1*grain.thickness*grain.atmosphere_drag_coeff_vert)* - abs(.5*atmosphere_curl - grain.ang_vel)* - (.5*atmosphere_curl - grain.ang_vel) + pi * grain.areal_radius^4. * rho_a * + (grain.areal_radius / 5. * grain.atmosphere_drag_coeff_horiz + + .1 * grain.thickness * grain.atmosphere_drag_coeff_vert) * + abs(.5 * atmosphere_curl - grain.ang_vel) * + (.5 * atmosphere_curl - grain.ang_vel) nothing end diff --git a/src/grid.jl b/src/grid.jl @@ -472,10 +472,10 @@ function conformalQuadrilateralCoordinates(A::Vector{Float64}, b = (delta*beta - alpha*epsilon) - (kappa*dx - gamma*dy) c = alpha*dy - delta*dx if abs(a) > 0. - d = b^2./4. - a*c + d = b^2. / 4. - a*c if d >= 0. - yy1 = -(b/2. + sqrt(d))/a - yy2 = -(b/2. - sqrt(d))/a + yy1 = -(b / 2. + sqrt(d)) / a + yy2 = -(b / 2. - sqrt(d)) / a if abs(yy1 - .5) < abs(yy2 - .5) y_tilde = yy1 else @@ -534,6 +534,7 @@ function findEmptyPositionInGridCell(simulation::Simulation, nx, ny = size(grid.xh) + i_iter=0 for i_iter=1:n_iter overlap_found = false @@ -917,9 +918,9 @@ function fitGridToGrains!(simulation::Simulation, grid::Any; max_x += padding max_y += padding - const L::Vector{Float64} = [max_x - min_x, max_y - min_y] - const dx::Float64 = 2.*max_radius - const n = convert(Vector{Int}, floor.(L./dx)) + L::Vector{Float64} = [max_x - min_x, max_y - min_y] + dx::Float64 = 2. * max_radius + n = convert(Vector{Int}, floor.(L./dx)) if 0 in n || 1 in n println("L = $L") println("dx = $dx") diff --git a/src/interaction.jl b/src/interaction.jl @@ -92,9 +92,9 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int) simulation.grains[j].poissons_ratio) # Effective normal and tangential stiffness - k_n = E*A_ij/R_ij + k_n = E * A_ij/R_ij #k_t = k_n*ν # Kneib et al 2016 - k_t = k_n*2.*(1. - ν^2.)/((2. - ν)*(1. + ν)) # Obermayr et al 2011 + k_t = k_n * 2. * (1. - ν^2.) / ((2. - ν) * (1. + ν)) # Obermayr 2011 else # Micromechanical parameterization: k_n and k_t set explicitly k_n = harmonicMean(simulation.grains[i].contact_stiffness_normal, diff --git a/src/io.jl b/src/io.jl @@ -5,6 +5,7 @@ if typeof(Pkg.installed("JLD")) == VersionNumber import JLD hasJLD = true end +using Compat.DelimitedFiles ## IO functions @@ -173,7 +174,7 @@ function status(folder::String="."; id = replace(id, "./", "") id = replace(id, r".*/", "") time_s = @sprintf "%6.2fs" data[1] - time_h = @sprintf "%5.1fh" data[1]/(60.*60.) + time_h = @sprintf "%5.1fh" data[1]/(60. * 60.) percentage = @sprintf "%3.0f%%" data[2] lastfile = @sprintf "%5d" data[3] if data[2] < 99. diff --git a/src/ocean.jl b/src/ocean.jl @@ -102,7 +102,7 @@ function readOceanStateNetCDF(filename::String) u, v = interpolateOceanVelocitiesToCorners(u_staggered, v_staggered) time = convert(Vector{Float64}, - NetCDF.ncread(filename, "time")*24.*60.*60.) + NetCDF.ncread(filename, "time") .* 24. * 60. * 60.) h = convert(Array{Float64, 4}, NetCDF.ncread(filename, "h")) e = convert(Array{Float64, 4}, NetCDF.ncread(filename, "e")) @@ -334,10 +334,10 @@ function applyOceanVorticityToGrain!(grain::GrainCylindrical, draft = grain.thickness - freeboard # height of submerged thickness grain.torque += - pi*grain.areal_radius^4.*rho_o* - (grain.areal_radius/5.*grain.ocean_drag_coeff_horiz + - draft*grain.ocean_drag_coeff_vert)* - abs(.5*ocean_curl - grain.ang_vel)*(.5*ocean_curl - grain.ang_vel) + pi * grain.areal_radius^4. * rho_o * + (grain.areal_radius/5. * grain.ocean_drag_coeff_horiz + + draft * grain.ocean_drag_coeff_vert) * + abs(.5 * ocean_curl - grain.ang_vel) * (.5 * ocean_curl - grain.ang_vel) nothing end diff --git a/src/packing.jl b/src/packing.jl @@ -36,8 +36,8 @@ function regularPacking!(simulation::Simulation, r_rand = 0. pos = zeros(2) h = .5 # disc tickness - dx = r_max*2.*(1. + padding_factor) # cell size - dx_padding = r_max*2.*padding_factor + dx = r_max * 2. * (1. + padding_factor) # cell size + dx_padding = r_max * 2. * padding_factor srand(seed) for iy in 1:n[2] @@ -69,9 +69,9 @@ r_j). function generateNeighboringPoint(x_i::Vector, r_i::Real, r_j::Real, max_padding_factor::Real) - R = rand()*(r_i + r_j)*max_padding_factor + 2.*(r_i + r_j) - T = rand()*2.*pi - return [x_i[1] + R*sin(T), x_i[2] + R*cos(T)] + R = rand() * (r_i + r_j) * max_padding_factor + 2. * (r_i + r_j) + T = rand() * 2. * pi + return [x_i[1] + R * sin(T), x_i[2] + R * cos(T)] end export poissonDiscSampling diff --git a/src/temporal.jl b/src/temporal.jl @@ -112,8 +112,8 @@ function findLargestGrainStiffness(simulation::Simulation) if grain.youngs_modulus > 0. k_n = grain.youngs_modulus*grain.thickness # A = h*r - k_t = k_n*2.*(1. - grain.poissons_ratio^2.)/ - ((2. - grain.poissons_ratio)*(1. + grain.poissons_ratio)) + k_t = k_n * 2. * (1. - grain.poissons_ratio^2.)/ + ((2. - grain.poissons_ratio) * (1. + grain.poissons_ratio)) else k_n = grain.contact_stiffness_normal k_t = grain.contact_stiffness_tangential diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl @@ -105,11 +105,11 @@ function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical, grain.lin_pos += grain.lin_vel * simulation.time_step + 0.5 * grain.lin_acc * simulation.time_step^2. + - 1./6. * d_lin_acc_dt * simulation.time_step^3. + 1. / 6. * d_lin_acc_dt * simulation.time_step^3. grain.ang_pos += grain.ang_vel * simulation.time_step + 0.5 * grain.ang_acc * simulation.time_step^2. + - 1./6. * d_ang_acc_dt * simulation.time_step^3. + 1. / 6. * d_ang_acc_dt * simulation.time_step^3. grain.lin_vel += grain.lin_acc * simulation.time_step + 0.5 * d_lin_acc_dt * simulation.time_step^2. diff --git a/src/util.jl b/src/util.jl @@ -18,8 +18,9 @@ Returns one or more random numbers from a power-law probability distribution. max_val::Number = 1.) val = ((max_val^(distribution_power + 1.) - - min_val^(distribution_power + 1.))*rand(dims) + - min_val^(distribution_power + 1.)).^(1./(distribution_power + 1.)) + min_val^(distribution_power + 1.)) * rand(dims) .+ + min_val^(distribution_power + 1.)) .^ + (1. / (distribution_power + 1.)) if dims == 1 return val[1] @@ -38,6 +39,6 @@ function harmonicMean(a::Number, b::Number)::Number if a ≈ 0. && b ≈ 0 return 0. else - return 2.*a*b/(a + b) + return 2. * a * b / (a + b) end end diff --git a/test/netcdf.jl b/test/netcdf.jl @@ -10,7 +10,7 @@ info("#### $(basename(@__FILE__)) ####") info("Testing dimensions of content read from Baltic test case") ocean = Granular.readOceanNetCDF("Baltic/00010101.ocean_month.nc", "Baltic/ocean_hgrid.nc") -@test ocean.time/(24.*60.*60.) ≈ [.5, 1.5, 2.5, 3.5, 4.5] +@test ocean.time / (24. * 60. * 60.) ≈ [.5, 1.5, 2.5, 3.5, 4.5] @test size(ocean.xq) == (24, 15) @test size(ocean.yq) == (24, 15) @test size(ocean.xh) == (23, 14) @@ -27,7 +27,7 @@ u1, v1, h1, e1 = Granular.interpolateOceanState(ocean, ocean.time[1]) u2, v2, h2, e2 = Granular.interpolateOceanState(ocean, ocean.time[2]) @test_throws ErrorException Granular.interpolateOceanState(ocean, -1.) u1_5, v1_5, h1_5, e1_5 = Granular.interpolateOceanState(ocean, - ocean.time[1] + (ocean.time[2] - ocean.time[1])/2.) + ocean.time[1] + (ocean.time[2] - ocean.time[1]) / 2.) @test u1 ≈ ocean.u[:, :, :, 1] @test v1 ≈ ocean.v[:, :, :, 1] @test h1 ≈ ocean.h[:, :, :, 1] @@ -36,7 +36,7 @@ u1_5, v1_5, h1_5, e1_5 = Granular.interpolateOceanState(ocean, @test v2 ≈ ocean.v[:, :, :, 2] @test h2 ≈ ocean.h[:, :, :, 2] @test e2 ≈ ocean.e[:, :, :, 2] -@test u1_5 ≈ (ocean.u[:, :, :, 1] + ocean.u[:, :, :, 2])/2. -@test v1_5 ≈ (ocean.v[:, :, :, 1] + ocean.v[:, :, :, 2])/2. -@test h1_5 ≈ (ocean.h[:, :, :, 1] + ocean.h[:, :, :, 2])/2. -@test e1_5 ≈ (ocean.e[:, :, :, 1] + ocean.e[:, :, :, 2])/2. +@test u1_5 ≈ (ocean.u[:, :, :, 1] + ocean.u[:, :, :, 2]) / 2. +@test v1_5 ≈ (ocean.v[:, :, :, 1] + ocean.v[:, :, :, 2]) / 2. +@test h1_5 ≈ (ocean.h[:, :, :, 1] + ocean.h[:, :, :, 2]) / 2. +@test e1_5 ≈ (ocean.e[:, :, :, 1] + ocean.e[:, :, :, 2]) / 2.