collision-5floes-normal.jl (10374B)
1 #!/usr/bin/env julia 2 using LinearAlgebra 3 4 # Check for conservation of kinetic energy (=momentum) during a normal collision 5 # between two ice cylindrical grains 6 7 verbose=false 8 9 @info "# One ice floe fixed" 10 sim = Granular.createSimulation(id="test") 11 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 12 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 13 Granular.addGrainCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose) 14 Granular.addGrainCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose) 15 Granular.addGrainCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose) 16 sim.grains[1].lin_vel[1] = 0.1 17 sim.grains[2].fixed = true 18 sim.grains[3].fixed = true 19 sim.grains[4].fixed = true 20 sim.grains[5].fixed = true 21 22 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 23 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 24 25 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 26 # should become more correct 27 28 Granular.setTotalTime!(sim, 10.0) 29 sim_init = deepcopy(sim) 30 31 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 32 Granular.setTimeStep!(sim, epsilon=0.07) 33 tol = 0.2 34 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 35 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 36 37 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 38 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 39 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 40 @test E_kin_rot_init ≈ E_kin_rot_final 41 @test 0. < norm(sim.grains[1].lin_vel) 42 for i=2:5 43 @info "testing ice floe $i" 44 @test 0. ≈ norm(sim.grains[i].lin_vel) 45 end 46 47 48 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 49 sim = deepcopy(sim_init) 50 Granular.setTimeStep!(sim, epsilon=0.007) 51 tol = 0.02 52 @info "Relative tolerance: $(tol*100.)%" 53 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 54 55 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 56 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 57 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 58 @test E_kin_rot_init ≈ E_kin_rot_final 59 @test 0. < norm(sim.grains[1].lin_vel) 60 for i=2:5 61 @info "testing ice floe $i" 62 @test 0. ≈ norm(sim.grains[i].lin_vel) 63 end 64 65 66 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 67 sim = deepcopy(sim_init) 68 Granular.setTimeStep!(sim, epsilon=0.07) 69 tol = 0.01 70 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 71 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 72 verbose=verbose) 73 74 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 75 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 76 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 77 @test E_kin_rot_init ≈ E_kin_rot_final 78 @test 0. < norm(sim.grains[1].lin_vel) 79 for i=2:5 80 @info "testing ice floe $i" 81 @test 0. ≈ norm(sim.grains[i].lin_vel) 82 end 83 84 85 @info "# Ice floes free to move" 86 87 sim = Granular.createSimulation(id="test") 88 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 89 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 90 Granular.addGrainCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose) 91 Granular.addGrainCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose) 92 Granular.addGrainCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose) 93 sim.grains[1].lin_vel[1] = 0.1 94 95 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 96 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 97 98 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 99 # should become more correct 100 101 Granular.setTotalTime!(sim, 40.0) 102 sim_init = deepcopy(sim) 103 104 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 105 Granular.setTimeStep!(sim, epsilon=0.07) 106 tol = 0.2 107 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 108 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 109 110 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 111 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 112 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 113 @test E_kin_rot_init ≈ E_kin_rot_final 114 for i=1:5 115 @info "testing ice floe $i" 116 @test 0. < norm(sim.grains[i].lin_vel) 117 end 118 119 120 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 121 sim = deepcopy(sim_init) 122 Granular.setTimeStep!(sim, epsilon=0.007) 123 tol = 0.02 124 @info "Relative tolerance: $(tol*100.)%" 125 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 126 127 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 128 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 129 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 130 @test E_kin_rot_init ≈ E_kin_rot_final 131 for i=1:5 132 @info "testing ice floe $i" 133 @test 0. < norm(sim.grains[i].lin_vel) 134 end 135 136 137 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 138 sim = deepcopy(sim_init) 139 Granular.setTimeStep!(sim, epsilon=0.07) 140 tol = 0.01 141 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 142 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 143 verbose=verbose) 144 145 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 146 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 147 @test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol 148 @test E_kin_rot_init ≈ E_kin_rot_final 149 for i=1:5 150 @info "testing ice floe $i" 151 @test 0. < norm(sim.grains[i].lin_vel) 152 end 153 154 155 @info "# Adding contact-normal viscosity" 156 @info "# One ice floe fixed" 157 sim = Granular.createSimulation(id="test") 158 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 159 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 160 Granular.addGrainCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose) 161 Granular.addGrainCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose) 162 Granular.addGrainCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose) 163 sim.grains[1].lin_vel[1] = 0.1 164 sim.grains[1].contact_viscosity_normal = 1e4 165 sim.grains[2].contact_viscosity_normal = 1e4 166 sim.grains[3].contact_viscosity_normal = 1e4 167 sim.grains[4].contact_viscosity_normal = 1e4 168 sim.grains[5].contact_viscosity_normal = 1e4 169 sim.grains[2].fixed = true 170 sim.grains[3].fixed = true 171 sim.grains[4].fixed = true 172 sim.grains[5].fixed = true 173 174 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 175 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 176 177 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 178 # should become more correct 179 180 Granular.setTotalTime!(sim, 10.0) 181 sim_init = deepcopy(sim) 182 183 184 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 185 sim = deepcopy(sim_init) 186 Granular.setTimeStep!(sim, epsilon=0.007) 187 tol = 0.02 188 @info "Relative tolerance: $(tol*100.)%" 189 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 190 191 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 192 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 193 @test E_kin_lin_init > E_kin_lin_final 194 @test E_kin_rot_init ≈ E_kin_rot_final 195 @test 0. < norm(sim.grains[1].lin_vel) 196 for i=2:5 197 @info "testing ice floe $i" 198 @test 0. ≈ norm(sim.grains[i].lin_vel) 199 end 200 201 202 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 203 sim = deepcopy(sim_init) 204 Granular.setTimeStep!(sim, epsilon=0.07) 205 tol = 0.01 206 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 207 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 208 verbose=verbose) 209 210 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 211 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 212 @test E_kin_lin_init > E_kin_lin_final 213 @test E_kin_rot_init ≈ E_kin_rot_final 214 @test 0. < norm(sim.grains[1].lin_vel) 215 for i=2:5 216 @info "testing ice floe $i" 217 @test 0. ≈ norm(sim.grains[i].lin_vel) 218 end 219 220 221 @info "# Ice floes free to move" 222 223 sim = Granular.createSimulation(id="test") 224 Granular.addGrainCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose) 225 Granular.addGrainCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose) 226 Granular.addGrainCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose) 227 Granular.addGrainCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose) 228 Granular.addGrainCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose) 229 sim.grains[1].lin_vel[1] = 0.1 230 sim.grains[1].contact_viscosity_normal = 1e4 231 sim.grains[2].contact_viscosity_normal = 1e4 232 sim.grains[3].contact_viscosity_normal = 1e4 233 sim.grains[4].contact_viscosity_normal = 1e4 234 sim.grains[5].contact_viscosity_normal = 1e4 235 236 E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim) 237 E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim) 238 239 # With decreasing timestep (epsilon towards 0), the explicit integration scheme 240 # should become more correct 241 242 Granular.setTotalTime!(sim, 10.0) 243 sim_init = deepcopy(sim) 244 245 @info "Testing kinetic energy conservation with Two-term Taylor scheme" 246 sim = deepcopy(sim_init) 247 Granular.setTimeStep!(sim, epsilon=0.007) 248 tol = 0.02 249 @info "Relative tolerance: $(tol*100.)%" 250 Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose) 251 252 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 253 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 254 @test E_kin_lin_init > E_kin_lin_final 255 @test E_kin_rot_init ≈ E_kin_rot_final 256 for i=1:5 257 @info "testing ice floe $i" 258 @test 0. < norm(sim.grains[i].lin_vel) 259 end 260 261 262 @info "Testing kinetic energy conservation with Three-term Taylor scheme" 263 sim = deepcopy(sim_init) 264 Granular.setTimeStep!(sim, epsilon=0.07) 265 tol = 0.01 266 @info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)" 267 Granular.run!(sim, temporal_integration_method="Three-term Taylor", 268 verbose=verbose) 269 270 E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim) 271 E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim) 272 @test E_kin_lin_init > E_kin_lin_final 273 @test E_kin_rot_init ≈ E_kin_rot_final 274 for i=1:5 275 @info "testing ice floe $i" 276 @test 0. < norm(sim.grains[i].lin_vel) 277 end