tocean.jl - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tocean.jl (6438B)
---
1 #!/usr/bin/env julia
2
3 # Check if ocean-specific functions and grid operations are functioning
4 # correctly
5
6 @info "Testing regular grid generation"
7 sim = Granular.createSimulation()
8 sim.ocean = Granular.createRegularOceanGrid([6, 6, 6], [1., 1., 1.])
9 @test size(sim.ocean.xq) == (7, 7)
10 @test size(sim.ocean.yq) == (7, 7)
11 @test size(sim.ocean.xh) == (6, 6)
12 @test size(sim.ocean.yh) == (6, 6)
13 @test sim.ocean.xq[1, :, 1] ≈ zeros(7)
14 @test sim.ocean.xq[4, :, 1] ≈ .5 * ones(7)
15 @test sim.ocean.xq[end, :, 1] ≈ 1. * ones(7)
16 @test sim.ocean.yq[:, 1, 1] ≈ zeros(7)
17 @test sim.ocean.yq[:, 4, 1] ≈ .5 * ones(7)
18 @test sim.ocean.yq[:, end, 1] ≈ 1. * ones(7)
19 @test size(sim.ocean.u) == (7, 7, 6, 1)
20 @test size(sim.ocean.v) == (7, 7, 6, 1)
21 @test size(sim.ocean.h) == (7, 7, 6, 1)
22 @test size(sim.ocean.e) == (7, 7, 6, 1)
23 @test sim.ocean.u ≈ zeros(7, 7, 6, 1)
24 @test sim.ocean.v ≈ zeros(7, 7, 6, 1)
25 @test sim.ocean.h ≈ zeros(7, 7, 6, 1)
26 @test sim.ocean.e ≈ zeros(7, 7, 6, 1)
27
28 @info "Testing velocity drag interaction (static ocean)"
29 Granular.addGrainCylindrical!(sim, [.5, .5], .25, .1)
30 Granular.setTotalTime!(sim, 5.)
31 Granular.setTimeStep!(sim)
32 sim_init = deepcopy(sim)
33 sim.grains[1].lin_vel[1] = 0.1
34 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
35 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
36 Granular.run!(sim, verbose=false)
37 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
38 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
39 @test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
40 @test E_kin_lin_init > E_kin_lin_final # linear velocity lost due to ocean drag
41 @test sim.grains[1].ocean_stress[1] < 0.
42 @test sim.grains[1].ocean_stress[2] ≈ 0.
43
44 @info "Testing velocity drag interaction (static ice floe)"
45 sim = deepcopy(sim_init)
46 sim.ocean.v[:, :, 1, 1] .= 0.1
47 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
48 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
49 Granular.run!(sim, verbose=false)
50 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
51 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
52 @test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
53 @test E_kin_lin_init < E_kin_lin_final # linear vel. gained due to ocean drag
54 @test sim.grains[1].ocean_stress[1] ≈ 0.
55 @test sim.grains[1].ocean_stress[2] > 0.
56
57 @info "Testing vortex interaction (static ocean)"
58 sim = deepcopy(sim_init)
59 sim.grains[1].ang_vel[3] = 0.1
60 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
61 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
62 Granular.run!(sim, verbose=false)
63 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
64 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
65 @test E_kin_rot_init > E_kin_rot_final # energy lost to ocean
66 @test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
67 @test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
68 @test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
69
70 @info "Testing vortex interaction (static ice floe)"
71 sim = deepcopy(sim_init)
72 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
73 sim.grains[1].lin_pos[1] = 0.5
74 sim.grains[1].lin_pos[2] = 0.5
75 sim.ocean.v[1, :, 1, 1] .= -0.1
76 sim.ocean.v[2, :, 1, 1] .= 0.1
77 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
78 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
79 Granular.run!(sim, verbose=false)
80 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
81 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
82 @test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
83 @test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
84 @test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
85 @test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
86
87 sim = deepcopy(sim_init)
88 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
89 sim.grains[1].lin_pos[1] = 0.5
90 sim.grains[1].lin_pos[2] = 0.5
91 sim.ocean.v[1, :, 1, 1] .= 0.1
92 sim.ocean.v[2, :, 1, 1] .= -0.1
93 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
94 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
95 Granular.run!(sim, verbose=false)
96 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
97 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
98 @test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
99 @test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
100 @test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
101 @test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
102
103 sim = deepcopy(sim_init)
104 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
105 sim.grains[1].lin_pos[1] = 0.5
106 sim.grains[1].lin_pos[2] = 0.5
107 sim.ocean.u[:, 1, 1, 1] .= -0.1
108 sim.ocean.u[:, 2, 1, 1] .= 0.1
109 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
110 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
111 Granular.run!(sim, verbose=false)
112 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
113 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
114 @test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
115 @test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
116 @test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
117 @test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
118
119 sim = deepcopy(sim_init)
120 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
121 sim.grains[1].lin_pos[1] = 0.5
122 sim.grains[1].lin_pos[2] = 0.5
123 sim.ocean.u[:, 1, 1, 1] .= 0.1
124 sim.ocean.u[:, 2, 1, 1] .= -0.1
125 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
126 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
127 Granular.run!(sim, verbose=false)
128 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
129 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
130 @test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
131 @test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
132 @test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
133 @test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained